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ABSTRACT 


Two theoretical methods and the development of a guidance, navigation, and 
control rapid protoyping system address the issue of considering the integral partici- 
pation of feedback early in the design process. The first method addresses the problem 
of sizing the horizontal tail on a statically unstable transport aircraft. Dynamic con-’ 
straints, which implicitly involve the flight control system, may prove more restrictive 
than traditional static constraints. Recovery from a severe angle of attack excursion, 
or penetration of a vertical wind-shear, are formulated in terms of the solution to 
a convex minimization problem utilizing LMIs, and numerically solved. The second 
method addresses the problem of tracking inertial trajectories, with applications for 
unmanned air vehicles. This problem is posed and solved within the framework of 
gain-scheduled control theory. This leads to a new technique for integrated guidance 
and control systems, with guaranteed performance and robustness properties. Finally, 
a rapid prototyping system for the flight testing of guidance, navigation, and control 
algorithms for unmanned air vehicles is designed. ‘The system affords a small team the 
ability to take a new concept in guidance, navigation, and control from initial concep- 
tion to flight test. A proof-of-concept demonstration of the system is detailed when 
the new integrated guidance and control algorithm previously described is tested in 


flight on an unmanned air vehicle. 
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ie OVERVIEW 


A. LAYOUT OF THE CHAPTERS 


Methods tn integrated plant-controller design and integrated guidance-controller 
design are contained in three chapters and two appendices. Applicable supporting 
computer code is explained in the appendices. The organization and content of the 
body of the work is as follows. 

Chapter II addresses the problem of sizing the horizontal tail on a statically 
unstable transport aircraft. A brief background in linear matrix inequalities prepares 
the reader for the problem formulation. Recovery from a severe angle-of-attack ex- 
cursion, or penetration of a vertical wind-shear, is formulated in terms of the solution 
to a convex minimization problem, and numerically solved. The effects of a num- 
ber of parameters in the aircraft definition, such as canards and flexible motion, are 
addressed. 

In Chapter III, the problem of unmanned air vehicles tracking inertial trajec- 
tories is addressed. The concept of trimming trajectories is defined, and the tracking 
of a trimming trajectory is converted into the problem of driving a generalized er- 
ror vector - which implicitly includes the distance to the trajectory - to zero. The 
linearization of the generalized error dynamics along the trajectory is shown to be 
time-invariant. Using these results, the problem of trajectory tracking is posed and 
solved within the framework of gain-scheduled control theory. This leads to a new 
technique for integrated guidance and control system design for unmanned air vehi- 
Glee 

The development of a rapid prototyping system for flight testing of guidance, 
navigation, and control algorithms for unmanned air vehicles is presented in Chapter 
IV. The system affords a small team the ability to take a new concept in guidance, 
navigation, and control from simulation to flight test. A proof-of-concept flight test 


demonstration of a new integrated guidance and control algorithm is detailed. The 


project includes a parameter identification problem for an unmanned air vehicle at the 


Naval Postgraduate School. Controller synthesis and implementation follow, based 


on the identified model. Finally, the system is used to test the theory presented in 
Chapter III in flight. 


B. 


CONTRIBUTION 


This section states the original contribution of this work. 


The problem of sizing a horizontal control surface for a statically unstable 
transport aircraft is addressed. Closed-loop performance of the plant and 
controller in response to certain dynamic constraints is used to define a design 
region that captures the integrated nature of the plant definition and controller 
synthesis process. 


A numerical tool is developed that determines acceptable combinations of 
tail volume, center-of-gravity location, and peak actuator rate for a statically 
unstable transport aircraft using state-feedback, or output-feedback, linear 
control to recover from a severe angle-of-attack excursion, or vertical wind- 
shear penetration. The method formulates the dynamic constraint as a convex 
optimization problem, and uses recently developed LMI solvers to solve the 
problem numerically. 


Numerical results demonstrate how to quantify the effect of the following in 
terms of their influence on the sizing of the horizontal control surfaces. 


— The addition of a canard to an aircraft definition. 
— The use of dynamic output versus static, state-feedback linear control. 


— The effect of simple, symmetric aeroelastic motion. 


Work on implementing nonlinear gain-scheduled controllers has been extended 
to trajectory tracking control problems for unmanned air vehicles. A broad 
class of trimming trajectories is defined, and a technique is developed that 
provides independent control of the space and time coordinates along trajec- 
tories in the class. Using this technique, it is shown that the performance and 
robustness of the linear design is recovered locally by the nonlinear plant and 
gain-scheduled controller. 


A unified system for the synthesis, simulation, and flight testing of avionics 
algorithms on unmanned air vehicles is constructed. The system is used to 
take new theory in integrated guidance and control from conception to flight 
test. 


bo 


ele INTEGRATED PLANT CONTROLLER 
OPTIMIZATION 


A. INTRODUCTION 


The use of an automatic flight augmentation system is commonplace on a 
modern aircraft. The benefits of its use may include such things as the remedy of 
undesirable flight characteristics, the reduction of pilot workload, and an increase in 
performance and fuel efhciency. Whether required for safe flight or not, it 1s hard 
to imagine a new transport aircraft being built without a sophisticated flight aug- 
mentation system. In fact, the next generation of high-speed civil transport (HSCT) 
aircraft will require some form of feedback control for safety of flight considerations, 
since current designs have an unstable short period. 

Since HSCT will require some form of automatic flight control always being 
active, the sizing of its control surfaces is no simple matter. Traditionally, static con- 
straints have been used to size the horizontal tail. For a given tail volume, constraints 
are calculated that limit the fore and aft travel of the center of gravity. Constraints 
that limit the forward center-of-gravity position include (1) sufficient nose-up pitch 
acceleration at the rotation speed (nose-wheel lift off), and (2) sufficient nose-up 
pitch acceleration at the approach speed in the landing configuration (go-around). 
Constraints that limit aft center-of-gravity position include (1) at brake release with 
maximum thrust, sufficient weight on the nose gear (tip back), (2) pitch-up acceler- 
ation at the rotation speed (nose-wheel lift off), and (3) sufficient nose-down pitch 
acceleration at minimum flying speeds [Ref. 26]. 

At the aft center-of-gravity locations projected for the approach flight condi- 
tion of the HSCT, dynamic constraints may be more restrictive than static constraints. 
The need to include dynamic considerations in the configuration definition process 
has been addressed before. References [Ref. 5] and |Ref. 4] describe early published 


work by Beaufrere in this area, largely motivated by the X-29 research program. In 


[Ref. 41], Schmidt uses the system sensitivity function to describe the fundamental 
trade-off that exists between the level of static instability that can be controlled and 
vehicle flexibility. Previous (unpublished) work in industry has used time domain 
analysis to determine how far aft (how unstable) the center of gravity could be before 
the airplane was unrecoverable. The analysis included a given design angle-of-attack 
disturbance, and given rate and position limits of the pitch control effector. 

This work extends the work of [Ref. 26], and uses a numerical technique similar 
to previous work in [Ref. 33]. There an integrated aircraft controller design method- 
ology using LMIs was applied to the control power sizing for an F—14 aircraft. The 
major contribution of this work is two fold. First, the tail sizing design problem is 
defined in terms that include the integral participation of a feedback control system. 
Since the degree of control of the longitudinal dynamics depends on how fast and far 
the longitudinal control surface(s) can be moved by the control actuator(s), a natural 
metric that captures the szze of the automatic flight control system is the maximum 
actuator rate. The design trade-off naturally includes consideration of actuator per- 
formance. For instance, it may be more cost effective to incorporate faster, generally 
larger and more expensive, actuators rather than pay the te penalty associated 


with a larger horizontal tail. In that light, we introduce the following definition. 
Tail Sizing Design Space: 


The Tail Sizing Design Space is the region of “acceptable” combinations of tail 
volume (Viz), center-of-gravity station (x29.), and peak actuator rate (timaz). 
The triplet {Vz 2.9.5 Umar}, defines an aircraft model and an automatic flight 
control system. The model is obtained through the linearization of the nonlin- 
ear dynamics of the HSCT at an equilibrium point. It is partially defined by 
the specified tail volume and center-of-gravity position. The automatic flight 
control system is characterized by the specified maximum actuator rate in the 
triplet. By “acceptable” it is meant that, for the model associated with the 
triplet Virece os ee a linear controller is known to exist that 1) stabilizes 
the plant, 2) meets prescribed dynamic performance constraints, and 3) does 
not exceed the maximum actuator rate in response to the dynamic constraint. 


Second, a numerical tool is developed that determines the Tail Sizing Design 
Space for a given HSCT dynamic model, flight condition, and dynamic constraint. 
This tool is termed the Tail Sizing Design Tool, and provides the capability to measure 
the effect of adding a second horizontal control surface in the form of a canard. The 
Tail Sizing Design Tool also provides the capability to measure the effect of simple, 
symmetric, flexible motion of the vehicle. Based on numerical analysis of designs, 
conclusions are drawn as to the relative effective of the use of canards, the use of 


static versus dynamic controllers, and the inclusion of aeroelastic effects. 


1. Example of a Design Space 

As a motivating example, consider a simple case. Let a(s) and b(h) be two in- 
dependent variables that are smooth functions of s and h, respectively. Furthermore, 
assume that they are independent of time, 1.e. Hels) = )and ath)) = 7 et 


d 

= = a(s)z(t) + b(h)u(t), (11.1) 
where z,u € WR’. Then, equation II.1 describes a family of dynamic systems. 
Evaluating (11.1) about arbitrary values of the parameters a(s) and b(h) results in 


the LTI systems 


&- 
o— 
oo 
— 
| 


a(so)x(t) + b(ho)u(t), 
=| 0) ern (11.2) 


It is desired to use feedback of the form 


to stabilize (11.2). Lyapunov stability theory states that this is equivalent to the 


existence of a p > 0 such that 





2p{a(so) — b(ho)k(so,ho)} <0. (1.4) 
This implies that 
k(so,ho) > ni (11.5) 


guarantees local stability of (II.1) at (so, ho). This places a lower bound on k(5o, ho). 
Assume that the absolute value of u(t) is constrained to be less than some 
nominal value, umaz- For this simple system, the maximum value of |u(t)| occurs at 


t = 0 and is equal to 
Umar = [A (so, ho)xo|. (11.6) 
This places an upper bound on k(S9, ho), 


k(so, io es (II.7) 


Lo 


The two limits, (II.5) and (11.6), are shown graphically in Figure 1. 


Region of 
Acceptable Gains 


os k(s,h) 


Figure 1. Range of Feedback Gains 


Uu 


If the boundaries, cat Sa amace allowed to approach each other, the region 


of acceptable feedback gains approaches a single point where 


a( 80) _ Umax 
nooo (11.8) 





Suppose a(s) varies linearly from a minimal value of 0.1 to a maximum value of 1.0. 
Suppose b(h) varies linearly from a minimal value of 0.3 to a maximum value of 0.9. A 
value of 1 was used for x9. The family of LTI systems, which validate the relationship 
in expression II.8, is represented by the curved surface shown in Figure 2. 

Below the surface, the LTI systems are not stabilizable subject to the stated 
constraint. Above the surface, multiple controllers exist that stabilize a given system 


subject to the stated constraint. If the maximum acceptable value of u(t) is included 
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Figure 2. Family of LTI systems with a single acceptable gain. 


as a horizontal plane, the Design Space for this simple triplet, {a, ,b , Umar}, consists 


of the region above the surface in Figure 2 and below this horizontal plane. This is 


shown in Figure 3. 





0 a(s) 


Figure 3. Example of a Design Space 


The intersection of the two surfaces is a line which represents limiting combi- 
nations of a(s) and 6(s) for a fixed Umaz and Zo. 

The chapter begins with a background in Linear Matrix Inequalities sufficient 
to formulate the general problem in a form suitable for solution numerically. Next, 
the Tail Sizing Design Tool is developed for a rigid-body HSCT model. The dynamic 
constraint associated with the problem is a severe angle-of-attack disturbance. Nu- 
merical results follow based on an application of the Tail Sizing Design Tool to an 
HSCT model representative of current designs. Then, a method is developed that 
models the effects of simple, symmetric, flexible motion of the HSCT. Numerical re- 
sults follow based on applying the Tazl Sizing Design Tool to an aeroelastic HSCT 
model. Comparisons to the rigid-body model are made. Finally, the Tazl Sizing De- 
sign Tool is developed for a rigid-body HSCT model where the dynamic constraint 
is the penetration of a vertical wind-shear. The chapter ends with conclusions and 


recommendations. 


B. BACKGROUND IN LINEAR MATRIX INEQUALI- 
TIES 


It is usually not possible to express the solution to a complex problem with 
multiple constraints analytically. Progress in computational power, however, make 
numerical solutions to these problems attractive. Efficient convex optimization algo- 
rithins, recently made available from The Math Works (Ref. 17], allow problems that 
can be formulated in terms of LMIs to be solved exactly. 

Two major results in LMI theory are drawn on to formulate the problem of 
sizing a longitudinal control surface. The first concept is that of an invariant ellipsoid. 


Consider an LT'I system: 


ee 


®’ 
| 
QQ 

= 


(11.9) 


where z € RR”, z € R?. Let P > 0 and define 
eae 0. Fe = |). 


Then € is called an invariant ellipsoid associated with (11.9) if, for every trajectory x 
satisfying (11.9), r(0) € € implies x(t) is in €, for allt >0 [Ref. 6]. 
Let € be an invariant ellipsoid for (II.9), and define V(r) = «7 Px. From the 
definition of €, it follows that V(x) < V(2(0)) < 1, and 
V(x) = 2 (F™P+PF)r <0, 
Se ee, 


On the other hand, suppose 3P > 0 such that F’P + PF <0 and 2x(0)? Px(0) <1. 
Then, for z(t) satisfying (11.9), we obtain 


x? (F7P+ PF)r =: V(r) <0. 
Integrating both sides from 0 to T' > 0 we get 
V(x(T))—-V(x(0)) < 0, 
= V(2(T)) < V(2(0)) <1, 
for any 7’ > 0. Therefore, we have shown that € is an invariant ellipsoid associated 
with the linear system (II.9) if and only if F7 P+ PF <0. 
Next, we will show how the idea of invariant ellipsoids can be. used to obtain 


bounds on the peak of the output in response to a known initial condition. Let € be 


an invariant ellipsoid associated with the linear system (II.9). Then, 
Olea eet) = max ('G'G¢ Wie A) 
Now, 


max ("TGC Serine pm (GaGa! e, 


[\u|]<1 
oa GP? |/?, 
Saath Guar. |): (11.10) 


where Amax(.) denotes the maximum eigenvalue of the argument. Therefore, an upper 
bound on the peak of |/z(¢)|| in response to x(0) can be found as a square root of the 


minimum of 6 subject to: 


$= \nae PCG Gian (11.11) 
2(0)’ Px(0) <1, (11.12) 
F'’P+ PF <0. (11 


Nonlinear inequalities, such as equation II.27, can be converted to LMI form using 


Schur complements [Ref. 22]. In general, any nonlinear inequality of the form 
R= 0,-1.0 = 6 Rescue ae (11.14) 


for symmetric Q and R, is equivalent to the LMI 


Qs 


0. (11.15) 
coe 








To see this, consider the congruence transformation, 


Ca) 


(11.16) 
SERIE: —R-'Ssi ] 0 R 


ie She 
0 Ee 


I ‘: : i 0 

















where the block diagonal structure makes equivalence to expression IJ.14 obvious. 


Using Schur complements on expression II.27 yields, 


&=\mnar( P98 4@ Ge) 


$1= PG GPa 


























6b Ghee I 0 O5eG ie) eea0 
= 0 << 
Pac eGe eee 07 Po | Gera 
6 G EG? 
2 > 
nO G 6 











In summary, the square root of a solution to the following optimization problem: 
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min 6 


subject to 
ieGae 
ean (1.17) 
G 6 
2(0)? Pz(0) < 1, (11.18) 
P’P+PRF < 0, (11.19) 
0): (11.20) 


provides an upper bound on the peak of ||z(t)|| for the system (II.9) in response to 
the initial condition x(0). 

The second concept used is that of an LMI region. In their recent work, 
Gahinet and Chilali [Ref. 7] introduce the concept of an LMI region, and have shown 
that such regions can be characterized by LMIs. LMI regions include sectors, disks, 
conic sectors, strips, as well as the intersection of any of these. For our purposes, the 
intersection of just two regions will suffice. The first region is a conic sector centered 
at the origin with half inner angle ¢. The second region is the left half-plane with its 
right boundary at —(@. The intersection of these regions is shown as the region CL in 


Figure 4. 


Im(s) 


Figure 4. Region L 


tl 


The eigenvalues of the system described by equation II.9 lie in the region £, 


if and only if there exists a positive definite matrix P, such that 


sind¢(F'P+ PF) cos¢(—PF+ F'P) 0 
cos¢(PF—F'P) sing(F’ P+ PF) 0 <0. (11.21) 
0 0 F'P+ PF + P2p 


1. State Feedback Synthesis LMIs 


A feedback interconnection of an LTI system and state-feedback controller is 


¢ =Ar+ BU 
zt =(A+ BK]z 
2 = eer 
i ine 
wherex € R", wE R™, z € KR. Substitution into the previous invariant ellipsoid 


and pole placement region matrix inequalities result in: 


Tinee ce 


subject to 


ie (C+ DK)? 


2), 
7 2 
(C+ DK) ae 
(ON) (Oe 
(sin ¢)(A’ P+ K? BT P4 PA+ PBK) (cos¢)(—PA— PBK +A! P 4K! BT p) 
(cos¢)(PA+ PBK —AlP—KTBT Pp) (sing)(A’ P+ K' BI P+ PA+ PBK) 
0 0 
0 
T ToT p. 0, 
AT p+KTBT P4 PA+ PBK + P28 


P > 0. (ies 


The above matrix inequalities are no longer affine in the unknown parameters, 
since the state-feedback controller, A’, and positive definite matrix, P, appear as 
multiplicative terms. Since P is positive definite, let Y = P~'. Perform the following 


congruence transformations, and perform the substitution of variables, W = WY. 


2 




















jg (Cac \* 
0< 
(C+ DK) ee 
i 0 PVC PDI) Y 0 
eo I 
0 Lf eS Pee Diab a: 
i (eS Diny > 
0< 
(C+ DK)Y Bo 
g (CY + DW)? 
ze (11.24) 
(CY + DW) ee ay 


earl 
ook 
oo 


0 (sin g)(A7 P+ KT BI P+ PA+PBK) (cos¢)(—PA— PBK + ATP 4 KTBT p) 
0 
Y 


(cos¢)(PA+ PBK —ATP—KTBT py) (sing) ATP 4+ KT BT P4 PA+ PBK) 
0 0 


0 Y 0 0 
0 0 Y 0 
Al P+K! BY P+ PA+ PBK + P26 0 o Y 


(sin o( YA! + Ww? BT 4 AY + BW) (cos¢)(-AY — BW +4 YA! 4 WT BT) 
(cos¢(AY + BW-—YA? —wl Bl) (sing yal + wT Bl +4 ay + BW) 
0 0 


0 
0 Q. 
YAl +wl Bl 4 ay + BW + Y28 me 


(11.25) 


Using Schur complements: 


IA 


z(0)* Px(0) > 


Osa 0) 8 1 


0. (11.26) 


Once again, the expressions I].24, I1.25, [1.26 constitute an LMI in the un- 
known matrix variables, W and Y. Suppose 3 Y > 0 and W such that they validate 
expressions 11.24, II.25, 11.26. Then, K = WY7~! is the state-feedback controller that 


maintains ||z|| < zmaz in response to the initial condition x(0), Vt > 0. 
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2. Output Feedback Synthesis LMIs 


A feedback interconnection of an LTI system and strictly proper, output- 


feedback controller is 








Ge Aes 

ee | A BCE 

2 = Cape BG ee (11.27) 
zy, = Apx, + Bry Zz = Cp apap: 1, 

Ue age 


where zs € R",c, €R™, WER”, ce Ro yee eee ama 
where 7? = [x7 z}]. The initial condition is 7(0)? = ng = [26 x4]. 
The four analysis matrix inequalities based on using the concept of an invariant 


ellipsoid and pole placement region are gathered for convenience. 


Pee 
0, (11.28) 
G 61 
vee Ue (11.29) 
caro eee (11.30) 
(sind)(F7P+PF) (cos¢)(—PF + F'P) 0 
(cos¢)(PF—F?P) (sind)(F'P + PF) 0 < 0 
0 0 FPP + PF + P2p 
(11.31) 


If the closed-loop system matrices in expression II.27 were substituted into 
expressions IJ.28, [[.29, II.30, and II.31, the expressions would not be affine in the 
unknown matrix vaniables. Similar to the state-feedback case, the key is to find a 
congruence transformation with a “linearizing” effect on the unknown matrix vari- 


ables. In a recent work by Scherer and Gahinet [Ref. 8], the authors prove that such 
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a transformation exists for the output feedback case. Consider the following partition 


eee and PP. 



































ee ay) 
a (11.32) 
eso 
7 R M 
et 3) 
M’ V 
Since P and P~! are both symmertic, so are R and S. From the identity 
I 0 
POE = 
0 7] 
- R M De Ay. 
Mi v N? Q 
RS+MN? RN+MQ 
= ; (11.34) 
WMIOS+VN?' MIN+VQ 
we can infer 
RS+MN* =1, 
or 
NM? =1-SR. (1E.35) 
Also, 
ees) I 
P = | (11.36) 
N? 0 
and 
R jf 
es = . (11.37) 
Mt 0 
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These identities are used to define the key matrices used in the congruence transfor- 


mation as follows: 







































































rey ol es: 
M’ 0 0 
=; Xo, (11.38) 
and 
Hs! maa 
Ce = : 
DEIN! MT 0 
= (11.39) 
Notice that 
= 
me: iS, ais: 
XoX7 = 
OPIN] S| |e Ciara 
= P. (11.40) 
Using X, and X2, consider the following congruence transformation: 
xX? 0 ||P GP eau MPG Sees 
| = (11.41) 
0 Tf G ol Carri GG 6] 
Expanding the (1,1) block of expression I1.41, we obtain 
KU PX, = 2G Coes Ge 
= ATX, 
RM a) 
ro {lo nTy 
ses (11.42) 
fa a 
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Expanding the (2,1) block of expression II.41 for the general case of output feedback, 


we obtain 








; cont 
GX = C. eC ) 
M? 0 
= Cri eC. CC, | (11.43) 
The matrix inequality in [1.28 becomes 
R fein Deen Ct 
I 5 ee = (0) (11.44) 


GR D,C.Mt G Ze 


mar 


After a congruence transformation with X,, the LMI in (II.29) becomes 


ie 
i 


EO: (11.45) 








This LMI is accounted for as a principle sub-matrix of the LMI in (1/.44). Schur 


complements are used to express the LMI in (II.30) as 


1 fa 
ore ro (1.46) 
no Po 


Assume zt, = 0, and consider the following congruence transformation: 


eG i 
no Po? 


0 x? 


1 0 
0 Xo 


! ng X2 
X? no De aXe 




















1 1g Xo 
Xi no DEG 


(11.47) 
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Expanding the (1,2) block of equation II.47 yields, 


i 
Ope 


T 


rf = | at eG | (11.48) 








After the appropriate congruence transformation, and substituting (11.48) into (II.47), 


we have 
eS 
Lo R I za (11.49) 
Seon I 5 


To see how the pole placement region LMI is transformed, look at the trans- 


formation of the (1,1) block of equation IT.31. 


FTP+PF < 0 <= 
POpffsztRPP < 0 <= 


Xd PEE Xx eee ae (11.50) 


Since P~! = X, Xz! and P-T = Xz" X7, equation 11.50 becomes 


XI Xo Xe (11.51) 


For output feedback synthesis 


A BC, 
B.C Ax 


F = (11.52) 








Substituting (11.52) into the first term in (II.51) results in 


I O 
met. 


Al eEGe 
BG een 


foe 
M* 0 


eG, 


} 
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A BC, R oI 
SA+NB.C SBC.+NA. || MT 0 | 
AR+ BC,MT A 


SAR+NB,CR+SBC,M' + NA,M!’ SA+NB,C 








(11.53) 


Define the change of variables, 


C, = C,M’, (11.54) 
eC. (11.55) 
A, = SAR+NB,CR+SBC,M' + NA,M', (11.56) 


and equation II.53 1s rewritten as 


=e 0) 


AR+ BC, A 
A; SA + B, 





Each term in the pole placement LMI] is transformed in a similar fashion. This 


results in 
sin @(=? +=) cos ¢(—H+ 57) 0 
cos@(={—=") sind(=’ +2 0 
( of < 0. (11.58) 
i 
0 0 SeP4+e+ 28 
I ] 








Notice that expressions I1.49, I1.45, 11.44, and IJ.58 are affine in the unknown 
matrix variables A;, B,, C,, R, and S. Thus, they constitute an LMI. Given 
Ar, Br, Cr, R, and S that validate expressions II.49, II.45, 11.44, and 11.58, the 
output-feedback controller is reconstructed as follows. First, reconstruct the matrices 


M and N via a singular value decomposition. 


HO 


Us aa ae (11.59) 


N-) = SV, (11.60) 
Ve (11.61) 


Form the controller by inverting the change of variables in expressions I1.55, I1.56, 


and I1.56 as 


A, = N"(Ay—SA,R—B.CpR—-— SB,Cy)M 7, (11.62) 
By =e (11.63) 
C, = C.M-, (II.64) 


HIGH ANGLE OF ATTACK RECOVERY 
1. Problem Formulation for a Rigid Body HSCT 


The Plant Controller Optimization (PCO) problem to be solved in this section 


can be stated as follows: 


Let an HSCT flight condition be specified by the aircraft’s flight speed, altitude, 
and flight path angle. Furthermore, let a certain set of flying quality require- 
ments exist for the HSCT at that flight condition. Let the dynamic constraint 
be defined as the recovery of the aircraft from a high angle-of-attack excursion, 
while not exceeding certain peak actuator rate and actuator amplitude limits. 
Then, for a range of horizontal tail volumes, and for a range of peak actua- 
tor rates, determine the aft center-of-gravity limits for which there still exists 
a linear feedback controller that 1) stabilizes the plant, 2) satisfies the flying 
quality requirements, 3) satisfies the dynamic constraint. 


Let xz., denote the center-of-gravity location as a fraction of the reference 


chord, and let x denote the vector of HSCT longitudinal states. Let u denote the 


horizontal tail incidence angle, and let Vy denote the tail volume. Figure 5 provides 


a description of the tail volume, Vy. 


This representation is convenient since, for one flight condition, the wing-body 


mean aerodynamic center remains constant. Thus, the control surface volume and 
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Figure 5. Definition of Vi. 


center of gravity are decoupled in terms of their influence on the vehicle’s dynamics. 
The longitudinal dynamics of the aircraft can be expressed by the following system 


of differential equations: 


ca eal ? ? Cc a, ? 
pre 7 Pw Mh) (11.65) 
z = Ett |e 


where F(.) is a C! function of x, u, t-9, Vx, H is a C! function relating x and z, and 
z = [Vr 7]? defines the true airspeed (Vr) and flight path angle (7). Let zo, ug denote 
the trim values of z and u for a given 2 and center-of-gravity location, 2.,,, 1.€. 
F (x0, Uo, Zego) Vp) = 0 and H(xo, uo) = 2. Here the vector z is used to characterize 
the flight condition. Later, when z,, is allowed to change, the trim values, ro and uo, 
will be recomputed for given values of x.,,, 29 and Vy,. Now, linearizing G about 


Lo, Ug yields a system of linear differential equations, 


6¢ = A(z, Feg,, Vi, Ot + B(Z0, Leqy, Vy OU, (11.66) 


Za 


where éz and éu denote small perturbations in z and u about Zp and wp, respectively. 
In the numerical analysis, the velocity components of 6x were normalized by the trim 
value of the true airspeed. 

Let the actuator dynamics be independent of (Zo, 2cg,, Vy, ), and let them be 


given by the set of differential equations, 


Oly = AOL au 
Gu ="Cxo7 (11.67) 
6 C07 ae ol, 
where éu denotes the actuator amplitude and éu denotes the actuator rate. Append- 
ing the actuator dynamics in series with the linearized longitudinal dynamics results 


in the system: 


A PC. 0 


a nee on 
: 0 Ay B, 
Gr (2020 VI=) 5 Tyo V6, (11.68) 
64 =| p Cmloveamiom 
thee 
fy = | 527 gat |. (11.69) 


Flying quality requirements are typically characterized by the level of attention 
and skill required of the pilot to control the aircraft. They are grouped in three levels. 
A lower level corresponds to more benign flight characteristics. In order to achieve 
certain Level II flying qualities requirements, the eigenvalues of G; must be placed 
in a more restrictive region in the left half plane. Figure 6 shows suggested locations 
of the Category B, closed-loop pole locations of the HSCT short-perod mode. These 
locations meet Level II flying quality requirements for the flight condition used in this 
study. Notice, these locations can be characterized as having a minimum damping 


ratio of 0.2 and natural frequency of 0.2 radians per second. This suggests that, in 
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order to meet Level II requirements for the HSCT, the short-period eigenvalues of 


the closed-loop system G; must be placed in the region in the complex plane given 


by Figure /. 


imaginary 





Figure 6. Acceptable Short Period Pole Locations 


Re(s) 





Figure 7. Region £ for Level I] Flying Qualities 
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2. Canard Configuration 

The horizontal tail is the most significant control effector for longitudinal con- 
trol of the aircraft. It is not uncommon, however, for long, slender aircraft designs 
to include canards for supplemental longitudinal control. The XB-70 and B-1B are 
two examples. Their presence is usually attributed to flying quality control issues 
involving the flexible nature of these aircraft [Ref. 1]. This point is addressed in a 
subsequent section concerning control of an aeroelastic aerodynamic model. To lay 
the ground work for that section, and as a baseline point of comparison for a rigid 
structure, this section addresses the addition of a longitudinal control effector for- 
ward of the wing. The effect on the Taz Sizing Design Space is explored when the 
feedback control system is free to utilize both control surfaces, in order to recover 
from angle-of-attack excursions. 

The revisions to the nonlinear equations of motion required to account for 
the addition of the canard parallel the development in [Ref. 14]. When these non- 
linear equations of motion are linearized about the equilibrium point determined by 


(20, Lcqy1 Vo Vey), the following LTI system results, 


Of = Aza. aaer Vito, Voy) 62 + Be( 205 Leqq, Voy )SUc + Br(Z0, Lego, Vp )OUn, (11.70) 


where u, and u; represent movement of the canard and horizontal tail, respectively, 
and Vg, is the canard volume. Actuator dynamics for both the canard and horizontal 
tail are appended to the linearized longitudinal dynamics as before. The actuator 


dynamics are given by: 
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Ola 


Mae 
du. 


dup 
dte 


dup 








This results in the following system, 


Gi (oe e553 Vito Ven) a 


where 








A,. 90 ies hema) dU- 
= OL, : 
0 Ae 0 Bar buh, 
1p 
= OTe 
== Oe aire eas 
Herr Ota 
= Og + DO, 
= Cone + 1D OTE. 
A a ee 
ieee dU. 
Oa 0 ou + 
0 sae Ou), 
0 ae 
dU; =|0 om ste 
dU. 
bu. =|0 Cx 0 |v +| D,. 0 
oun, 
iu, =[0 0 Ca, |v 
dU, 
sin =[0 0c, ]iv+[o v,, | | 
dup, 
(11.72) 
Te 
Oa | 527 oc 27, | 


The initial condition is defined by the angle-of-attack perturbation. Let the 


angle-of-attack perturbation be given as ao; then the initial condition is 


dV 


[0 V,cos~'(ag) 0 0 0 


a 


For each design point, G (20, Lego; Vito, Vey); the question becomes, is the set of feed- 


back controllers that recover the aircraft from the angle-of-attack excursion, while 


20 


maintaining acceptable flying qualities, and without saturating the actuator, or ex- 


ceeding a certain actuator rate, empty? 


3. Proposed Numerical Solution 

In this section, we make the conjecture that the problem at hand is analogous 
to the example presented in the introductory section. The center-of-gravity location 
plays the role of the parameter influencing the stability of the open-loop plant. The 
tail volume plays the role of the parameter influencing how controllable the plant is 
from the input where feedback is to be applied. It is shown that when the plant 
parameters are fixed, the problem can be formulated in terms of an LMI, which is 
affine in both the controller parameters and actuator limits. Minimizing the actuator 
rate limit is a convex optimization problem that can be solved exactly. For one set of 
plant parameters, let the result of the process be characterized by the final values of 
all of the parameters. Repeat the process for different sets of values of the plant pa- 
rameters. Then, an assumption 1s made that the points lie on a smooth hypersurface. 
If the grid of plant parameters is sufficiently fine, the shape of the hypersurface can 
be discerned. The projection of the hypersurface into the three dimensional space 
spanned by tail volume, center-of-gravity location, and actuator rate limit defines a 
surface in the Tail Sizing Design Space that is extremely useful in the design of the 
aircraft. 


Let 


, (G; (Zo, L ego 4 Va Veo ), ae Une rea 3 Uewae Vo) = 


Vee > 0% 
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0. C,, OY + [D,, OW 2, a 
l T 
(oo ag 
Vo \ 


Y tren iy (hp 





(00 C,,] ¥ + [0 D,,]W uy 
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where 
ee O.g ee LC a, 
2 _ i ee) 
A=]0 Ag. 0 = 
eee 
0 0 rien 
Let 


05(G, (Zo, Lego Vito 3 Vo, Ne lina ) eee Cat 3 i ae Vo) ae 


{ Ax, (SPC. 1h 
peer mG eC? 
(ers > 0, 
Co en. 2 
ee CeCe. 
aoe} > 0, 
C..Ck 0 wh... 


2d 


(11.73) 


(11.74) 


R I CRUNCH REC! 

I S ATICT 0] > 0, 
[0 C,.Aa] R+C,.BaCr [0 C,.] Ac Ue 

R I CPBICT + R'[ATC. 0] 

I S ATICT 0) > 0, 
(0 C,.Aa] R+C;.BaCk [0 C;,] Aa Wh mar 


[eran goes 


Zo R i = a 


nq TTT tS 


ee) I o 
(sin ¢)(117 + 1) (cos ¢)(—I1 + Il’) 0 
(cos¢)(I1I—II7) (sin ¢)(II? + IN) 0 ae 
RI _ 
0 0 7 +H+ 26 
eS 
(11.75) 
where 
A BeCae oe anean ea: AOBeGo. Pe aRean 
Wit. =e ec ae Jes | “ Bap | ox aM re 
aa ; | A BcCag BprCa, ; 
Ar S| 0 Aac 0 + By, 
0 0 Aa; 
(11.76) 


Then, a sufficient condition for the existence of a dynamic, output feedbac- 
controller, or static, state-feedback controller, that stabilizes the feedback system 
Gi (Zo, X ego VHo, Vey), does not exceed actuator amplitude and rate limits, and results 
in acceptable flying qualities in response to the angle-of-attack excursion defined by 
Yo, is for the sets, ®,, or ®2 to be non-empty. At this point, we observed that the 
same numerical results were obtained using a slightly different approach. Instead of 
minimnzing a parameter Umaz in the LMI decision vector, we exploited the fact that 
the short period pole becomes more unstable as the center of gravity is moved aft. 


This allowed the use of a feasibilty algorithm vice minimization algorithm to acheive 
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the same results. The benefit was a significant reduction in computational time. Now, 
the Plant Controller Optimization (PCO) problem considered in this section can be 


stated as follows: 


Maximize {z,,} 

Subject to: 

F (Xo, Vo, Leg, Vz, Vo.) = 9, 

H(X0, Uo) = Zo, 

(Y,W) € ®1(G: (Leg, VGA nal, eatin lee 3 te, Vo), ei 
Or 

(Az, Bu, Cr, R,S) € ®2(Gi (teg, Vito» Veo); Umar > Uhmarr Ucmar> Ucmas YO): 


(11.78) 


A solution to this PCO problem includes a linear controller that stabilizes the feedback 
system G (Zo, 2c, VHp, Ve,), and meets actuator limit requirements, as well as a 
providing a maximum aft center-of-gravity location. 

The numerical solution used to map the Tail Sizing Design Space involved a 


binary search over center-of-gravity stations: 


1. Fix Vi, Ve, Veet ee — land z= 0. 

re hg et Deg )/ 2. 

Ol Ce lexanc Gtesn)— Uy Mt(.Xq, U5) = Zo for Xo, Up: 

4. Obtain A( Xo, Vo, Leg,), B(X0, U0, Leg.) and form G; (Zo, Leg, )- 

5. Solve for (Y,W) € ®1(G; (2eq), Uhmars Uhmar> Ucmar> emar> Yo) OF; 

6. Solve for (Ax, By, Cy, R, S) € ®2(Gt (cq), Umass Uhmaz> Ucmar > teemass YO): 
7. If no such (Y, W) exist 


e vcgmaz a Lego) 


else 


Zo 


~ vegmin = Lego: 
8. If fegngs — Lean Ors BOO 


9. Increment Vy, go to 1. 


4. Results - Rigid Body HSCT 

This design methodology was applied to a number of aerodynamic models 
representative of current high-speed civil transport designs. The results are shown 
for the aerodynamic model termed Ref A (see Appendix A). Appendix A details how 
a convenient build-up of the nonlinear longitudinal dynamics in terms of wing-body 
and tail contributions facilitates the process of mapping the Tail Sizing Design Space. 

The state feedback synthesis LMIs (set ®,) were used. Figure 8 shows an 
optimization run for a single tail volume of 0.2, no canard, and peak actuator rate 
limit of 30 degrees per second. As the center-of-gravity location is moved aft, the peak 
actuator rate approaches the limit. The actuator amplitude remained well below its 
saturation limit, which was set at 20 degrees of travel from trim. 
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Figure 8. Fixed Tail Volume and Peak Actuator Rate Limit 


Figure 9 details the process for one tail volume and a sweep of peak actuator 


rates. Again, the tail volume was 0.2 and the peak actuator rate limit was incremented 
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from 5 to 30 degrees per second. Notice, initially the center of gravity moves quickly 
aft with only small increases in the peak actuator rate required. However, at some 
point, the peak actuator rate required becomes very sensitive to changes in the center- 
of-gravity station. For the Ref A model, this break point is in the vicinity of 0.565 


percent mean aerodynamic cord. 
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Figure 9. Fixed Tail Volume, Sweep of Peak Actuator Rates 


This process is repeated for the range of tail volumes of interest. Figure 10 
shows the results of this process for a tail volume of 0.1, and a tail volume of 0.2. 
A two dimensional representation of the data 1s cumbersome. When repeated for 
numerous tail volumes, the Tazl Sizing Design Space could be viewed in three dimen- 
sions. Figure 11 shows the result of fitting a surface to data obtained for a range of 
tail volumes from 0.1 to 0.3, and a range of peak actuator rates from 5 to 30 degrees 
per second. The surface represents a lower bound on the peak actuator rate required 
by the feedback control system to recover from the angle-of-attack excursion, for var- 
ious combinations of center-of-gravity location and tail volume. An upper bound in 
the Tail Sizing Design Space would be a fixed limit on the peak actuator rate avail- 


able. Figure 12 shows the inclusion of this plane for a value of 20 degrees per second. 
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The volume below the plane and above the curved surface represents the Tail Sizing 
Design Space, where linear state-feedback controllers are known to exist that meet 
design requirements. Figure 13 shows the intersection of the plane representing the 
peak actuator rate available of 15 degrees per second, with the curved surface rep- 
resenting the minimum peak actuator rate required for a rigid-body HSCT with no 
canard. This is seen to be an aft center of gravity limit line on the standard center 


of gravity versus tail volume (scissors) plot. 
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Figure 10. Two Tail Volumes, Sweep of Peak Actuator Rates 


At this point, the canard volume was fixed at 0.05. Following the procedure 
outlined earlier, the lower surface of the Tazl Sizing Design Space was mapped for 
a range of horizontal tail volumes from 0.1 to 0.3, and for a range of peak actuator 
rates from 5 to 40 degrees per second. Amplitude and actuator rate limits for the 
canard and horizontal tail were matched. Figure 14 shows the resulting lower bound 
in the design space. Figure 15 shows the results for Ref A HSCT with and without 
the canard. Figure 16 shows a slice of the two surfaces in Figure 15 at a peak 
actuator rate limit of 15 degrees per second. This should be familiar as a conventional 


scissors plot. Of course, as one would expect, the design space increases with the 
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additional control power available from the canard. The pertinent question is whether 


or not there is any benefit in spreading the control power for and aft, so to speak, 


or if the Tail Sizing Design Space would have increased just as much had the tail 
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volume alone been increased by 0.05, and the canard not added. In general, it makes 
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Figure 12. 3D Tail Sizing Design Space with Peak Actuator Rate Limit 
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Figure 13. 2D Slice of the Tail Sizing Design Space 


more sense to compare the competing configurations in terms of equal amounts of 
combined horizontal tail and canard surface area. For instance, profile drag is more 
closely related to the surface area of the control surfaces, among other things, and 
the design goal might be to minimize drag for the same aft center-of-gravity station 
and actuator rate limit. For this example, which utilized the Ref A data, the distance 
from the vehicle’s wing-body neutral point to the aerodynamic center of the canard 
or horizontal tail was the same. Therefore, comparisons in terms of normalized area 
or volume are equivalent. 

Figure 17 compares the two configurations, the first without a canard and the 
second with a canard. The percent change in total control volume required in going 
from a configuration without a canard to a configuration with a canard is shown as a 
function of the aft center-of-gravity limit. The same flying quality requirements and 
actuator limits were used. For the rigid-body HSCT model, the benefit gained from 
the inclusion of a canard is about a 10 percent savings in total control volume. 

This process was repeated, using the output feedback synthesis LMIs (set ®2). 


Figure 18 shows the resulting Taz Sizing Design Space obtained utilizing dynamic, 
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the plane representing a target rate limit of 25 degrees per second is shown in Fig- 
ure 20. By inspection, it 1s clear that the limits of the Design Space are nearly 


full-state, feedback controllers. Figure 19 compares these results with those obtained 
utilizing static, state-feedback controllers. The intersection of the two surfaces with 
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Figure 15. Tail Sizing Design Space With and Without a Canard 
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Figure 16. 2D Slice of Tail Sizing Design Space With and Without a Canard 


identical. This is an important and somewhat surprising result considering the added 


complexity of the output-feedback controller. 
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Figure 17. Decrease In Total Control Volume Through the Addition of a Canard 
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Figure 18. 3D Tail Sizing Design Space - Dynamic Controller 


Problem Formulation - Aeroelastic Model 


The development of an integrated aeroelastic aerodynamic model is well docu- 


mented in work by Waszak and Schmidt [Ref. 45]. The nonlinear equations of motion 


were derived using a Lagrangian approach and some standard simplifying assump- 
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Figure 19. Tail Sizing Design Space Comparison: Static versus Dynamic Controllers 


3f 


tail volume 





“10 20 30 50 60 70 


40 
c.g. (%c) 
Figure 20. 2D Comparsion: Static versus Dynamic Controller 


tions. The elastic deformation is assumed sufficiently small such that linear elastic 
theory holds. Informally, the flexible deformation is captured by the contribution of 
an infinite number of states termed the generalized elastic coordinates. Associated 
with each generalized elastic coordinate is an zn vacuo mode shape. In laymens terms, 
this describes the shape of a single elastic mode when fully deflected. Of course, for 
practical purposes, only a finite number of elastic modes are retained for analysis. In 
a body-fixed reference coordinate system, the resulting elastic airplane equations of 


motion for the longitudinal dynamics are: 


Qs, 
Q:, 

Lyyq = Qo, 
(hy Gert) ee ae 


m(u + qw + gsiné) 


m(w — qu — gcos@) (11.79) 


where, 
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Q; = generalized forces, 
m; = generalized mass, (11.80) 
ni; = generalized elastic coordinates. 


Let X, and Z be the total aerodynamic and propulsive forces along each of the 
axis in the body-fixed reference frame. Similarily, let M be the total aerodynamic 
and propulsive moment about the body-fixed y-axis. Define 6z, éz to be virtual 
displacements along the z and z axis, and define 6@ be a virtual rotation about 
the y axis. When combined with the virtual pertubations of the elastic generalized 
coordinates, 6n;, the virtual displacements are called the generalized coordinates. The 
Principle of Virtual Work is used to expand the generalized force terms in (II.80) 
through the use of the generalized coordinates. It states that 
O(déW) 


Q; = Aba)” (11.81) 


where 6W is the work associated with an arbitrary virtual displacement of the gen- 





eralized coordinates, 6q;. This virtual work done by the aerodynamic and propulsive 
forces, relative to the inertial reference frame, expressed in a coordinate system at- 


tached to the body of the aircraft, is: 


éW = Xbx+ 262+ [M +(2X —2Z)] 0p + [ P(x,z)-S>$;6n:dS. (1.82) 
ill 


The last term represents the work done by the distributed surface pressure 
due to virtual displacements of all of the elastic generalized coordinates. 
Applying (II.82) to the expression for virtual work, the generalized forces are 


seen to be 








Q.= a(6z) = X, (11.83) 
— O(6W) | 
Q.= A(z) > Z, (11.84) 





OC — = M +(zX —2xZ), (11.85) 


—aswy a 2 
Qn = Fay = Hea I P(x, z)- 6:67:48). (11.86) 


At this point, a method needs to be chosen to determine the aerodynamic and 





propulsive forces and moments. We assume that a knowledge of wing-body, tail, and 
canard aerodynamic stability and control derivatives for the rigid-body aircraft are 
available. Furthermore, the vehicle is assumed to have one flexible mode, with its 
mode shape and generalized modal mass known. Generally, the first symmetric mode 
shapes of all the HSCT designs are very similar. Each mode increases the number of 
states in the linear model by two, which subsequently increases computational times 
required by the Taal Sizing Design Tool. Since the modes retained in the linear model 
are those that we wish to actively control, the maximum number of modes retained 
would be three or four. This method can be combined with other methods, [Ref. 45], 
and used to capture the residuals from the truncated modes. The method used here 
works well for capturing the interaction of the rigid-body states with the flexible- 
body states for the first few symmetric modes. It is worth noting that the LMI based 
Tail Sizing Design Tool tool can utilize any method that is capable of generating a 
linear aeroelastic model at a specified center of gravity, tail volume/canard volume, 
and flight condition. The main contribution of this method is its suitability to the 
iterative nature of the numerical solution, and its use of widely available rigid-body 
aerodynamic stability and control derivatives. 

Therefore, the aerodynamic and propulsive forces are expressed in terms of a 


body-fixed reference frame as 


X = Lsina— Dcosa+ 7;, 
Z = —Lcosa— Dsina+ Ty, (11.87) 


where L and D are the lift and drag aerodynamic forces, a is the angle of attack, and 


T; is the component of thrust in the z,, direction. The total lift of the airplane is 


AQ) 


ee Lay, be Le (11.88) 


where the subscripts are meant to be suggestive of the canard, wing-body, and tail 


contributions. In coefficient form, the expression becomes, 


aoe 


CL 5 


S 
Cmte Ore (11.89) 


where 5S, S;, and S, are the reference area of the wing, tail and canard respectively. 
For the problem at hand, (II.89) is rewritten explicitly in terms of a tail and canard 
volume. Let c denote the reference mean aerodynamic cord, J; denote the distance 
from the aerodynamic center of the tail to the aerodynamic center of the wing-body, 
and |, be a similar length associated with the canard. Then, (II.89) can be expressed 


in terms of control volumes as 


Cc. = VorCu, + Om. + Vasu. (11.90) 
c t 


The individual lift coefficients are expanded in a first order Taylor series ex- 
pansion about the components of the state vector and control inputs. Local changes 
in the free stream dynamic pressure have been incorporated into the local lift-curve 


slope derivatives. For example, the coefficient of lift for the tail is expressed as follows: 


u ; 
of Chi, + Cir. V + Cy,.a+ Ci, & + C14 
V 
SECU + 2 Chie ete 2 Chin, ni + oon Lye (11.91) 
Most of the terms in (II.91) are familiar as aerodynamic stability derivatives 
for a rigid-body aircraft. Some, such as Cy,,, may be new to the reader. In order to 


derive an approximate expression for these aeroelastic stability derivatives, we need 


41 


to introduce the concept of a mode shape. It is assumed that the general elastic defor- 
mation of the unconstrained body can be described as the product of two functions, 
one a function of only spatial coordinates, and the other a time dependent function. 
Informally, the mode shapes, or free vibration modes, are the spatial function, and 
the generalized elastic coordinates are the time function. Then, the local relative elas- 
tic displacement is described in terms of the mode shape, ¢;(xr), and the generalized 


coordinate, 7;(t), as 


l= gia (11.92) 





7 = = Ee) ni(t). (11.93) 


Expressions [1.92 and II.93 can be used to obtain formulae for the force and 
moment dependence on flexible motion. Again, considering the tail contribution to 
lift as an example, the lift-curve slope dependence at the tail on 7 is computed from 
the rigid-body aerodynamic stability and control derivatives, and from expressions 


II.92 and II.93 as 





ie 
Cie = ere 7 oe (11.94) 
Cry, = Cr, a. (195) 


where Zac, is the normalized location of the mean aerodynamic center of the tail. The 
canard and wing-body terms are treated in an analogous fashion. Once the total lift 
is calculated, an assumed knowledge of a drag polar is used to calculate the total 
drag. 

Similar to the build-up of the total lift on the vehicle, the total pitching mo- 
ment on the vehicle is expressed in coefficient form. Again, a first order Taylor series 


expansion about perturbations of the state variables and contro] inputs results in 
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Cu = Cup + Cu.o + Cue + Cuyg + Cay, Ur 


+Cu.te + >) Ca,, 1 + >) CM, Hi. (11.96) 


—|) t=1 
The individual derivatives are written in terms of the center-of-gravity location 
(x.,), aerodynamic center of the wing-body (2ac,,), and familiar control volume terms 


[Ref. 14]. For instance, 

















Cm. = Cx1,(2eg — Lacy,) — Vit a 
Lh se ql mite (11.97) 
and 
Oye Cr (i — lace.) — Vit ae 
+Ve = 2 - =i (11.98) 


which are composed of stability and control derivatives that were either assumed 
provided or previously computed. 
The remaining term to consider involves expansion of the generalized force 


associated with the elastic generalized coordinates. Recall, 


e : 
Qn = FEylf, Ple=): bidndS). (11.99) 


As a first step, the integral expression in (11.99) is separated into three parts: 


i" P(z,2)-¢,6n,dS = | P(x, z)- ;6n.dS + . P(a,z)- di6nid5 
i / P(x, z)- $;6n;dS. (11.100) 
t 
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The pressure distribution over each surface 1s approximated by a point force acting at 
the aerodynamic center of the lifting surfaces of the aircraft, and the axial component 


is ignored. Therefore, the three integrals above are approximated as 


[ Pe. z) : b,dn:dS mS Fz, bi{ Lac. ON: aia Fz,,,0i( Lacy, ON: ae F'z,6i( Lac )ON- 


(Toa 
Then, using (I1.99), a reasonable approximation of the generalized forces is 
Ory ~ Fz,¢:( Zac.) aig Fz,,,0i(Zacys) a5 F7z,0i( ac,)- (11.102) 


Each term in (II.102) is expanded in a first order Taylor series expansion. The force 


on the tail, for instance, becomes 


U . 
7, = faa G + a as 7, Oe eee aU 
V 


a SS Zan Ni + yy FT, i + Zu, (i. 


2=1 t=] 


(11.103) 


Recall, the generalized forces, Q,, drive the dynamics of the elastic generalized coor- 


dinates. The rigid-body states couple into the elastic states via terms like 


PZ ad Oh aes ences aqSVin C1, ; (Lac, )a, (11.104) 
t 
and the control inputs couple into the elastic states through terms like 
Fa, Oi(taa)Ua = cos agSVin Cy, bi( Tacx Un: (11.105) 
t 


The remaining terms couple the elastic states among themselves. Using the 


fifth term in (11.103) for an example, we obtain 
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oS Veae, Pi( Dace), (11.106) 





Lae = COS aqS Vin CL., ( 
: t 


aud 


1Z,, = Cos aqSVin=Cn, §i(2ae,) eee 
: t 


(11.107) 

This completes the description of the modeling of the nonlinear aeroelastic 
dynamics. The intent was not to document each term. That level of detail is left to 
Appendix A. Instead, the flavor of how to incorporate the influence of a few elastic 
modes was demonstrated. The method assumed a knowledge of wing-body, tail and 
canard aerodynamic stability and control derivatives for the rigid-body, and a set of 
mode shapes and generalized modal masses. 

The nonlinear aeroelastic equations of motion were linearized at an equilibrium 


point determined by the flight condition, canard volume (possibly zero), tail volume, 


and center of gravity. The resulting linear model is given by: 


a A A as B Rie 
: _ | Arr Are ie eal | (11.108) 
LE Agr AEE LE BE Uh 
where 
T 
ZR = | w gq 0 : 
r 


Generalized elastic coordinates are somewhat lacking in physical intuition, and 
are not directly measured by any sensor suite. The relationships between sensed pitch 
angle and pitch rate at a given local body station (subscript /), the rigid-body states 


(subscript &), and the generalized elastic coordinates are: 


oO dd; 
b= Oe — (am (11.109) 
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er (11.110) 


=o; 
= _— —),,1;- Il. 
in Gt (11.111) 
Let 
ney 0422 
T=|0010 oOo (#), Ue 


define a similarity transformation that replaces the generalized elastic coordinates, 
7, 7, With a more physically meaningful pair of states, such as q and 6;. The location 
chosen for the placement of the local pitch rate and pitch angle sensors was the cockpit 
station. Typically, an area is sought where the mode shape slope has its most positive 


value {Ref. 1]. The linearized aeroelastic model now becomes 


Sr re es B 
Uc 
ir _ | ARR ARE | “ LT R | 
Agr AEE BE Uh 
0; 0; 
LR 
Uc 
=e | a | | | (11.113) 
Un 
0; 


Using Ref A data from Appendix A, the behavior of a rigid-body HSCT aircraft 
is compared with an aeroelastic model with the first elastic mode retained. Table I 
compares the eigenvalues of the two models. Figure 21 highlights the undesirable 
effect of the elastic motion. It shows the pitch rate sensed at the cockpit to a step 
input of the elevator for the two models. 

The frequency separation between the elastic mode and the short period dy- 
namics is approximately 1 Hertz, and the damping of the flexible mode is only 0.02. 
Typically, attempts are made to attenuate the feedback prior to excitation of the 


flexible dynamics. On large transport aircraft with the flexible dynamics close in 
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[_[_Short Period [Long Period [Flexible Mode 
[| frequency [ damping [frequency | damping | frequency [ damping 
[Aeroelastic] 0.90 | 75 | ol4 | 055 | 67 | 0.02 
[Rigit-body| 0.91 | 75 | O14 [007 | nja [aja 


Table I. Eigenvalues of Aeroelastic and Rigid-body Models 
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Figure 21. Pitch Rate at Cockpit, Aeroelastic vs. Rigid Body 


frequency to the short period dynamics, this is hardly possible. Furthermore, even 
if suitable notch or low-pass filtering within the control loop could be attained, the 
extremely light damping of the flexible modes results in problems eras of gust- 
induced structural responses and fatigue life. Therefore, the problem posed will be 
one of actively controlling the flexible modes retained. Generally, this will entail im- 
proving the damping of the flexible dynamics, while ensuring stability of the short 
period dynamics as the center of gravity is moved aft. 

As before, actuator dynamics are appended to the aeroelastic model. The 
Plant-Controller Optimization problem formulation is exactly the same as discussed 


in section l. 
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6. Results - Aeroelastic Model 

The results are shown for the aeroelastic dynamic model termed Ref B in 
Appendix A. Ref B utilizes the same rigid-body stability and control derivatives as 
Ref A. The effect of a single, symmetric, flexible mode was incorporated. The mode 
shape is shown in Figure 22. It was experimentally determined that requiring the 
feedback controller to increase the damping of the flexible mode to 0.05 resulted 
in acceptable damping of the short period mode. Once again, the LMI based Tail 
Sizing Design Tool was used to map out a surface in the Tail Sizing Design Space. 
Figure 23 shows the design space for Ref B without a canard. An upper bound in 
the Tail Sizing Design Space is shown by the level plane at a peak actuator rate of 
25 degrees per second. Figure 24 compares the aeroelastic model to the rigid-body 
model. Obviously, considerably more actuator rate is required for the aeroelastic 
model in order to recover from the angle-of-attack excursion, while actively controlling 
the flexible dynamics. Intuitively, this seems obvious. The Tail Sizing Design Tool 
provides a metric for comparison. For example, consider a target tail volume of 
0.20 and center-of-gravity station of 45 percent mean aerodynamic cord. This is well 
within the Tail Sizing Design Space of the rigid-body model, one outside that of the 
aeroelastic model. The aeroelastic model requires either an extra 0.05 increase in tail 
volume, or 14 degrees per second increase in peak actuator rate, to move within its 
Tail Sizing Design Space. In some respects, this is the price to be paid for actively 
controlling a non-rigid vehicle. 

Next, the effect of the addition of a canard to the aeroelastic model is explored. 
As before, a canard volume of 0.05 was selected, and the design space was determined. 
The rest of the design requirements remained unchanged, and the results are shown 
in Figure 25. Figure 26 compares the design space of the aeroelastic model with 
and without the canard. Figure 27 shows the resulting aft line on the tail volume 
sizing plot when the two are cut at a peak actuator rate of 25 degrees per second. 


Of course, the added control volume due to the canard allows a further aft center-of- 
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Figure 23. 3D Tail Sizing Design Space: Aeroelastic Model 
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Figure 24. Tail Sizing Design 
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Figure 26. Tail Sizing Design Space Comparison: Canard On and Off 


(Ref. 1]. This method provides a metric to quantify that benefit. In this example, 


the inclusion of a canard is 300 percent more effective when added to the aeroelastic 


model then when added to the rigid-body model. 
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As before, the process was repeated using the output feedback LMIs. 
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Figure 28. Decrease In Total Control Volume Through the Addition of a Canard 


results utilizing the aeroelastic model support the conjecture that dynamic controllers 
of the same order as the plant do no better than static, state-feedback controllers. 
A representative comparison is shown in Figure 29. There the aft center-of-gravity 
stations are shown for a range of tail volumes and fixed canard volume. The peak 


actuator rate limit was 40 degrees per second. 


D. GUST RECOVERY 


A second dynamic requirement imposed on transport aircraft is the ability 
to recover from a severe gust. The rationale is that the open loop HSCT should 
not be so unstable that an unrealistically fast actuator and flight control system are 
required. The gust recovery criterion is similar to that of the high angle-of-attack 
excursion in many respects. Both are at their most adverse when the vehicle is at a 
slow speed. Also, the aft center-of-gravity configuration is limiting. However, unlike 
the high angle-of-attack excursion, which is assumed to occur essentially instantly, 
gust recovery criterion involves the vehicle penetrating a particular wind-shear over 


time. 
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Figure 29. 2D Tail Sizing Design Space Comparison: Static versus Dynamic Con- 
troller on an Aeroelastic Model 


The wind-shear occurs in the vertical plane, begins at zero velocity, builds to 
some peak value, and then dies out again. Since the shear is defined to occur as 
the vehicle crosses into a fixed region in space, the rate at which it is penetrated 
depends on the flight speed of the vehicle. In some industry documents, the gust 
profile is defined to reach its peak within 30 feet of travel of the vehicle. However, 
there is ambiguity as to the rate and peak value of the gust. The bottom line is that 
relatively late in the design/testing process, a flight simulation will be conducted. 
The vehicle will penetrate some gust, based on empirical weather data gathered from 
airports throughout the country. Failure at this stage involves a major redesign. 
What is required is the ability to quantify the concepts used in the definition of the 
design requirement, such as so unstable, and unrealistically fast actuator. The Tail 
Sizing Design Tool provides a bench mark to assess competing configurations that is 


directly related to recovery from a gust penetration. 


1. Additional LMI Considerations 
Recall the state feedback synthesis LTI system description. 
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and suppose that the desired block structure of the feedback gain, KX, 1s 
ine Ky 0 | (11.115) 


A sufficient condition to assure the proper block structure of the feedback 
gain is to specify an appropriate block structure to the unknown variables in the 
corresponding LMIs. Recall, for the state synthesis feedback problem, the feedback 


gain is recovered as K = WY. Therefore, specify 
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2. Problem Formulation 

The gust recovery requirement describes the face of the gust profile as a One 
Minus Cosine disturbance. An acceptable alternative to the One Minus Cosine profile 
is to fit two exponential functions to the profile. This captures the important rapid 
change in airmass velocity as the vehicle penetrates the shear. Unlike the One Minus 
Cosine disturbance, the gust velocity returns to zero after reaching its peak. The 
rate of decay of the gust velocity after reaching its peak, however, can be made quite 


slow. Therefore, the impact on the vehicle’s dynamics is minimal. 
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Let t, be the time it takes the HSCT to penetrate the shear, and let Gp be 


the peak vertical velocity of the shear. Then, the following functions define the two 


wind-shear profiles. 


One Minus Cosine Profile: 


G f 
el cos(——_ )), Ve 
) t, 
= Gp, ee 


Double Decaying Exponential Profile: 
f(t) = —G, exp" +Gy exp", 0<t, 


where 


A 
A = ole and Az = AG: 


Both wind-shear profiles are shown in Figure 30. 





: 2 
t = distance/Vt 


Figure 30. Two Wind shear Profiles 
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Let the vehicle’s linearized longitudinal dynamics be given by the LT] system 


D0 


py = Av+ Bu, 


[| A, Ay Aq Ag | Aace |v + Bu, (11.121) 
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as before. Based on (11.119), the state space description of the wind-shear is 
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where v, is the vertical velocity of the shear. Now, a model of the vehicle’s dynamics 


in an airmass whose vertical velocity is v,(t) is 





A (A, sin 0) + Aw cos 99) (Aysin 89 + Ay cos Io) B 
7) — 1) = U, 
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The velocity states in (II.123) are inertially referenced quantities. Consider a 


similarity transformation based on the change of variables, 


¢ =fn, 


where 
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and 


¢ = TAT 1C+TBr, 
=: AC + Bu, 
Co = In. (11.123) 


Then, the first two states of ¢ are perturbations in true airspeed and angle of 
attack, both airmass relative quantities. This makes the first four states of ¢ natural 
quantities to be considered for feedback. However, for all practical purposes, the gust 
states are not measured, and should not be considered for feedback. This implies that 


the feedback considered should be of the form: 


= K, 0 | e | , (11.124) 


where (2 corresponds to the gust states, and ¢; corresponds to all of the remaining 


states. 


3. Proposed Numerical Solution 
Let 
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Now, replace ®, with 3 in (II.77), and proceed with the algorithm outlined 


in section (3). 


4. Gust Recovery Results 

This design tool was being developed concurrently with a NASA led multi- 
corporate effort to define a baseline configuration for an HSCT vehicle. The design 
engineers desired the ability to quickly asses the influence of changes in key parameters 
of the wind-shear on the final design. This necessitated that computational time be 
kept short, on the order an hour or less, since the intention was to investigate multiple 
scenarios throughout the day. 

In that light, the gust recovery criterion was applied to the numerical model, 
Ref A, with a fixed tail volume of 0.11 and no canard. Still to be finalized were the 
design requirements concerning the definition of the wind-shear. The design team 
desired to know how sensitive their designs were to changes in peak values of the 
wind-shear and/or time to penetrate the wind-shear. The Tail Sizing Design Tool 
was used to compute aft center-of-gravity stations for differing wind shear profiles. 

Figure 31 shows the sensitivity of the aft center of gravity location to changes 


in the peak value of the gust. The penetration distance is one cord length. The peak 
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gust velocities were 30, 45, and 60 feet per second. The data indicates that at around 
20 degrees per second peak actuator rate, every 50 percent increase in the peak gust 
velocity pushes the aft center-of-gravity station forward 10 percent. 
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Figure 31. Increasing the wind-shear peak velocity moves the center of gravity limit 
forward. 

On the other hand, figure 32 shows the sensitivity of the aft center of gravity 
location to changes in the penetration distance. The peak gust velocity was 45 feet 
per second. A penetration distance of 1, 13, and 2 cord lengths is shown. Choosing 
an actuator rate limit of 20 degrees per second again, every 50 percent increase in 
wind-shear penetration distance pushes the aft center-of-gravity station aft 5 percent. 

On a percentage basis, the aft center of gravity limit is about twice as sensitive 
to changes in the peak wind-shear velocity as it is to changes in the penetration 


distance. 


E. CONCLUSIONS 
The sizing of the horizontal control surface(s) for HSCT is a difficult problem. 
Traditionally, aircraft definition has taken place apart from feedback considerations. 


Controller design was a derivative process, with little contribution to the aircraft 
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Figure 32. Increasing the wind-shear penetration distance moves the center-of-gravity 
limit aft. 


definition. Of course, the aircraft were statically stable, and feedback control was 
viewed as a significant enhancement but not flight critical. This is an interesting 
example of an application where traditional methods simply lack the necessary tools 
to adequately define the aircraft configuration. Clearly, the aircraft definition will 
come first. However, it must be done with tools that provide a quantifiable knowledge 
of the impact on the demands of the feedback control system. 

The most intuitive metric to select, capturing the demands of the feedback 
control system, is the peak actuator rate required. The inclusion of this metric 
adds an extra dimension to the tail sizing problem. The two dimensional tail sizing 
scissors plot is inadequate, and a natural extension is the Tail Sizing Design Space. 
A conventional scissors plot can be recovered by viewing the intersection of the Tail 
Sizing Design Space with a level plane. However, the three dimensional space allows 
the designer to see how sensitive the aft boundary is to changes in actuator rate. 
Clearly, there is value in knowing when small changes in the center-of-gravity location 
result in large changes in the actuator rate required. 


While not solvable analytically, new efficient algorithms make the problem 
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tractable numerically. A new design tool is provided that provides the capability 
to quickly determine the Tad Sizing Design Space for a given aircraft configuration. 
This chapter demonstrates that many tail sizing problems can be formulated as LMIs, 
and exploit new interior point algorithms to obtain exact numerical solutions. The 
iterative nature of the design process requires a tool that can generate solutions in 
a timely manner. Using the Jail Sizing Design Tool, the user can make adjustments 
to the aircraft definition, and quantifiably asses the impact on the demands of the 
feedback control system. 

In the first half of the chapter, two fundamental changes in aircraft definition 
were explored. Their influence on the Jail Sizing Design Space for a representative 
model of an HSCT was quantfied. The first was the use of canards in addition to a 
horizontal tail. The second was the effect of simple, flexible motion. This resulted in 
the testing the matrix of basic aircraft configurations shown in table E. Additionally, 
the effect of using a static, state-feedback controllers, or dynamic, full-state, output- 


feedback controllers, was investigated for each element in the matrix. 


Rigid Body-Tail Only Rigid Body-Tail with Canard 





Flexible Body-Tail Only | Flexible Body-Tail with Canard 


Table II. Matrix of Aircraft Definitions Tested 


Numerical results suggest that canards provide a small benefit for a rigid body 
HSCT. Their use is more effective when used on a flexible body HSCT. The metric 
used to asses their effectiveness was the change in total horizontal control volume. 
In this example, an aeroelastic model realized a reduction in total horizontal control 
volume of approximately 30 percent through the use of a canard. The rigid-body 
model realized a reduction in total horizontal control volume of approximately 10 
percent. The use of a dynamic, full-state, output-feedback controller did not improve 


the results. 


6] 


In the second part of the chapter, a different dynamic constraint from the high 
angle-of-attack excursion was addressed. ‘The dynamic constraint was the penetration 
of a severe vertical wind-shear. The definition of the wind-shear profile was open. 
Therefore, the effect of changes in the wind-shear profile on the Tail Sizing Design 
Space for a rigid-body HSCT was explored. In this example, a boundary of the Tail 
Sizing Design Space was found to be twice as sensitive to changes in the peak level of 


the gust then to changes in the distance to penetrate the gust. 


62 


II. INTEGRATED GUIDANCE AND 
CONTROL 


A. INTRODUCTION 


In a great number of envisioned mission scenarios, Autonomous Vehicles will 
be required to follow inertial reference trajectories accurately in 3-D space. See, for 
example, [Ref. 36, 9, 11] and the references therein. Similar requirements emerge 
from the recent work at NASA on the descent trajectory synthesis for air traffic 
control [Ref. 43]. To achieve this goal, the following systems must be designed and 
implemented on board autonomous vehicles (AVs): 1) navigation, to provide estimates 
of linear and angular positions and velocities of the vehicle, 11) guidance, to process 
navigation/inertial reference trajectory data and output set-points for the vehicle’s 
(body) velocity and attitude, and ii1) control, to generate the actuator signals that 
are required to drive the actual velocity and attitude of the vehicle to the values 
commanded by the guidance scheme. 

The advent of GPS (Global Positioning System) has afforded AV systems en- 
gineers a powerful new means of obtaining accurate navigation data that is required 
for precise tracking of given inertial trajectories. However, traditional guidance and 
control schemes used to steer the vehicle along such trajectories may prove inade- 
quate in the case where frequent heading changes are required, or in the presence of 
shifting wind [Ref. 35]. Traditionally, such systems are designed separately, using 
well established design methods for control, and simple strategies such as line of sight 
(LOS) for guidance. See [Ref. 20] and (Ref. 29] for interesting applications to un- 
derwater and air vehicles, respectively. During the design phase, the control system 
is usually designed with sufhciently large bandwidth to track the commands that are 
expected from the guidance system. However, since the two systems are effectively 
coupled, stability and adequate performance of the combined system about nominal 


trajectories are not guaranteed [Ref. 35]. In practice, this problem can be resolved by 
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judicious choice of guidance law parameters (such as the ” visibility distance” in LOS 
strategy), based on extensive computer simulations. Even when stability is obtained, 
however, the resulting strategy leads to finite trajectory tracking errors, the magni- 
tude of which depends on the type of trajectory to be tracked (radius of curvature, 
vehicle’s desired speed, etc.) [Ref. 35]. 

This chapter proposes a new methodology for the design of guidance and con- 
trol systems for AVs, whereby the two systems are designed simultaneously. This 
methodology has two main advantages over traditional ones: 1) the resulting trajec- 
tory tracking system achieves zero steady state tracking error about any trimming 
trajectory, and ii) the design methodology explicitly addresses the problem of stability 
of the combined guidance and control systems. 

Earlier work based on this approach can be found in (Ref. 15, 42, 25]. In par- 
ticular, the authors introduce a methodology for the design of controllers for UAVs 
to track inertial trajectories that are given in space and time coordinates. The tra- 
jectories considered are equilibrium (also known as trimming) trajectories of AVs, 
which are helices parameterized by the vehicle’s linear speed, yaw rate and flight 
path angle. Furthermore, in [Ref. 15, 42, 25] it is shown the linearization of the vehi- 
cle error dynamics and kinematics about any trimming trajectory is time-invariant. 
Thus, the problem of designing integrated guidance and control systems for AVs to 
accurately track trimming trajectories can be solved by using tools that borrow from 
gain scheduling control theory, particularly those reported in [Ref. 27]. Within the 
framework of {Ref. 27], the vehicle’s linear speed, yaw rate, and flight path angle play 
the role of scheduling variables that interpolate the parameters of linear controllers 
designed for a finite number of representative trimming trajectories. The results in- 
troduced in [Ref. 27] on the D-implementation of gain scheduled controllers can then 
be used to obtain a combined guidance/control system such that the properties of the 
linear designs are recovered locally, about each trimming trajectory. An interesting 


and very important consequence of the D-implementation is that it leads naturally to 
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a controller structure where the only exogenous commands required are the desired 
linear inertial position and the yaw rate, thus avoiding the need to feedforward the 
trimming conditions for the remaining state variables and control inputs. 

However, the technique presented in [Ref. 15, 42, 25] has a shortcoming which 
may be of concern for UAVs and AUVs when tracking trimming trajectories in the 
presence of changing air and water mass. Since the controllers described in those 
references achieve accurate tracking of trajectories defined in terms of space and 
time coordinates, the relative speed of the vehicle with respect to the air cannot be 
controlled externally, its value being computed internally as a function of the tracking 
error. In practice, this may lead to unacceptable performance in the presence of winds, 
since the change in the airspeed of the vehicle may result in stall or structural damage. 

Clearly, eliminating the time coordinate in the trajectory definition and using 
the vehicle attitude to null out trajectory errors while maintaining constant airspeed 
should resolve this problem. A similar approach has been introduced in a number of 
publications on robot control. Of particular interest is the work reported in [Ref. 40}, 
where the subject of path following control for wheeled robots is addressed. See also 
(Ref. 35] for a detailed analysis of the stability of an autonomous underwater vehicle 
about nominal trajectories in the horizontal plane. 

In this chapter, these ideas are formalized within the basic framework for 3-D 
trajectory tracking controller system design developed in [Ref. 15, 42, 25]. Using 
the concepts outlined in [Ref. 40], the linear position of an AV is given in terms 
of its location with respect to the closest point on a desired trajectory, together 
with the arc length of an imaginary curve traced along that trajectory. Tracking of a 
trimming trajectory by the vehicle at a fixed speed is then converted into the problem 
of driving a generalized error vector - which implicitly includes the distance to the 
trajectory - to zero. Moreover, it is shown that the linearization of the generalized 
error dynamics about the corresponding trimming path 1s time-invariant. Using these 


results, the problem of trajectory tracking is posed and solved in the framework 
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of gain scheduled control theory, leading to a new technique for integrated design 
of guidance and control systems for AVs. The chapter summarizes the resulting 
design methodology, and illustrates its application to the design and implementation 
of a nonlinear trajectory tracking controller for the UAV Frog [Ref. 25]. Numerical 
simulations using a full set of nonlinear equations of motion of the vehicle show the 
effectiveness of the proposed techniques. 

The subject of trajectory tracking has also been addressed in the literature 
~ on control of nonholonomic vehicles. In [Ref. 44], the problem of tracking a nominal 
trajectory by a nonlinear system is considered. The key idea includes linearizing the 
nonlinear system along the trajectory, then using the resulting time varying lineariza- 
tion to obtain a time varying, state-feedback controller that locally exponentially sta- 
bilizes the system along the trajectory. The paper includes examples of applications 
of the proposed technique to a mobile robot and a front wheel drive car. The nominal 
trajectories considered in [Ref. 44] are not restricted to be trimming trajectories. 
However, all the examples presented consider the case of trimming trajectories only. 
Another approach is used in [Ref. 19], where a tracking problem for a surface ma- 
rine vessel is considered. Here the authors use feedback finearizaienieitn dynamic 
extension to obtain a controller to track trajectories that consist of lines and arcs of 
circles (a special case of trimming trajectories in the plane). 

The solution to the trajectory tracking problem proposed in this chapter dif- 
fers considerably from the ones introduced in [Ref. 44, 19]. Here, the key idea is to 
reduce the problem to the design of a tracking controller for a linear-time invariant 
plant utilizing a simple nonlinear transformation that inverts the vehicle kinematics. 
This poses no robustness concerns since kinematics are usually well known, partic- 
ularly in the case of air and underwater vehicles. It is important to point out that 
the application of the nonlinear transformation results in a nonlinear plant, whose 
linearization along trimming trajectories is time-invariant. Once in the linear setting, 


the designer is free to choose his favorite control synthesis technique to achieve the 
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desired closed-loop performance and robustness. The chapter provides a simple al- 
gorithm for implementing the linear controller on the nonlinear plant such that the 
properties of linear controller are preserved along each trajectory. This is in contrast 
to the approach in [Ref. 44], where the problem is reduced to that of designing an 
exponentially stabilizing state-feedback controller for a linear time-varying system. 
This leads to a controller design that is problem specific and does not address the 
issues of performance and robustness. On the other hand, in [Ref. 19] the authors 
point out that extending to air vehicles the feedback linearization technique use for 
trajectory tracking of the surface craft is difficult due to the unstable zero dynamics 
that are characteristic of aircraft. 

The chapter is organized as follows. Section B develops the rigid body dy- 
namics and introduces appropriate notation. Section C introduces the error dynam- 
ics necessary to solve the problem. Section D formulates and solves the problem of 
integrated guidance and control of UAVs for the class of trajectories introduced in 
Section C. Section E describes an application to the integrated guidance and contro] 


of the UAV Frog. Finally, Section F contains the main conclusions. 


B. RIGID BODY DYNAMICS 


This section begins with a review of the equations governing the motion of a 
rigid body, and specifically of an unmanned air vehicle. Notation familiar in robotics 
[Ref. 10], but not commonly used in the field of aircraft dynamics, is introduced. It 
is hoped that the familiar environment of the equations of motion of a rigid body will 
make the reader comfortable with the notation. It will be used extensively later in 
developing an integrated guidance and control algorithm. 

Let {J} denote an inertial reference frame. For the purposes of this devel- 
opment, the rotation of the earth and its associated Coriolis’ force can be ignored. 
Thus, a local tangent plane reference frame will be considered an inertial reference 


frame. Let {B} denote a body coordinate frame that is fixed with respect to the body 
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of the vehicle. Finally, let {W} denote a stability coordinate frame aligned with the 
velocity vector of the vehicle. 
Free vectors can be expressed in whatever reference frame is most convenient. 


The velocity of the vehicle with respect to {I} is a free vector. When resolved in 


{B}, it will be described by 
5 
vo= E v w | : ( Tes 


Free vectors can be transformed from one frame to another via the rotation matrix 
describing the relative orientation of the two frames. For example, it is common to 
describe the orientation of {B} with respect to {I}, by the Euler angles (¢, 0, w). 
Let 


then 7 R(A), or simply 7? R, is the rotation matrix from {I} to {B}, and given by 


B cos(@)cos{y) . cos(@) sin(w) . —sin(@) 
I R = sin(¢) sin(@) cos(w) — cos(¢) sin() sin(~) sin(@) sin(w) + cos(¢) cos(y) sin(¢) cos(theta) 4 
cos(¢@) sin(@)cos(y) + sin(¢)sin(y~) cos(¢) sin(@) sin( pst) — sin(¢) sin( pst) cos(¢) cos(@) 


(III.2) 


As an example, consider the position P of the vehicle. Its velocity, $P, in {7} can 


be expressed as 


Another common rotation matrix relates the orientation of { B} to the stability 
axis of the vehicle {W}. The vehicle’s angle of attack, a, and side-slip angle, £, 


define the orientation of {W} with respect to {B}. The associated rotation matrix 


is termed 7,R(a, 3). As an example, if V is the velocity of the vehicle with respect 
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to {I}, resolved in {B}, then 


WV || 
0 | = BR(a,B)V 
0 


is the same vector resolved in {W}. 

The angular velocity of {B}, with respect to {/}, resolved in {B} will be 
represented by 9. It is the set of familiar body-fixed rotation rates, roll rate (p), 
pitch rate (q), and yaw rate (r). The body-fixed rotation rates are related to the 


Euler angles by 


p I 0 — sin(0) = 
q| = | 0. cos(¢) — sin(¢) cos(6) “0 |, 
i 0 -—sin(¢) cos(¢) cos(@) £ 
= Q(A)A. (111.3) 


Inverting (III.3), the time rate of change of A is related to 9 by the matrix Q(A), or 


simply Q, as 


d 

a ele) 

a an’ 

iO I. (11.4) 


Using this notation, an application of Newton’s Law to the linear and angular 


motion of the vehicle yields 


d F 

hee) 7 

dt ao m’ 

d 

ae eee, ae NV, (ies) 


where F and WN are the total external force and moment acting on the vehicle resolved 
in {B}, mis the body mass, and Zz is the inertia tensor of the body resolved in {B}. 
Standard practice is to separate out the contribution of gravity, (G’), from the 


force and moment terms. The remaining terms are expanded in a first order Taylor 
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series expansion about some choice of state variables and control inputs. A convenient 


choice of state variables is (u, 8, a, p, q, r), and standard control inputs are elevator 


(de), aileron (6a), rudder (ér) and throttle (é6¢). Applying this to (III.5) results in 


d l 
—V = -Q.xV+4? RG + —pl|V||’sR.(a, 8) 
dt 2m 
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Cp, Cp, 0 Cp. U 
Cy eas 0 Cy, 0 B 
Cin Cy, © Gm Ox 
ét 
0 0 
de 
Cy es , (111.6) 
da 
0 0 
| or }J 
0 C1, 
0 ienlae (3) CMs 
sb Cn, 
b 
0 Ci, IMI 0 0 
Ca,0 0 am 9 | 
b 
0 Cn. 0 0 IVT 
; (II1.7) 


where 


S := wing reference area, 
C := wing mean chord, 
m =") ilass, 
Ibs := _aircraft’s inertia tensor resolved in {B}, 
p ‘= _ air density, 
b [= wing span, 
D, Y, L := — normalized drag, lateral force and lift, 
l, M, N := — normalized roll pitch and yaw moments, 
ce. := normalized nominal force or moment coefficient, 
Ox, = stability derivative:5*. 


Notice that some symbols in (III.6) and (III.7) have multiple meanings, and attention 
must be paid to the context in which they are used. The terminology is standard in 
aircraft dynamics. 

A convenient grouping of terms in (III.6) and (III.7) is to let H be a vector 
valued function of the control inputs, Z be a vector valued function of all the terms 
affine in #H, and F be a vector valued function of all the remaining terms. Then, 
equations III.6 and III.7 can be compactly expressed in terms of these functions. 
Subscripts are used to indicate the correspondence of the function to the appropriate 


linear or angular motion equation. As an example, 


l 
ly := 5 AllV I"sRu(a, 8) 0 0 Cy, CY, |- (111.8) 


Then, the complete equations of motion of the vehicle are expressed in this notation 


(6 


as 


2 V = Fy(V, Vad) ee cy eee 


40 = Fe(V,0,A) + Za(V,Q)H(V,,U), 


Go “ (111.9) 
dt ie =p rove 
f2£A =O%). 


dt 


C. GENERALIZED ERROR DYNAMICS 
1. Trimming Trajectories 


The objective of the guidance and control systems is to steer an autonomous 
vehicle along prescribed inertial trajectories Po € R°. Furthermore, we require that 
the vehicle be trimmed along any such trajectory. Let {C’} define the coordinate 
system attached to the vehicle, and let Ag define the desired inertial orientation of 
{C'}. The coordinate system {C’} represents the desired inertial orientation of the 
vehicle along Po. Therefore, at trim {B} = {C}. Next, we define the set € of the 


trimming trajectories about which the vehicle is expected to operate: 


Po £ Pe = RVo, 
Ea e 


= Ag = O(Ac) Ne =: Wee, 
Fy(Ve,%e, Xe) + iviVe, Me Ve, Va, 0) = 
Fa(Ve,Qe,Ac) + Za(Ve, 2Qc)H(Ve, 2c, U-) = 0, 





(IIT.10) 
where Vc, Qc and U, denote the trimming values of V,Q and U, respectively, and Ag 
denotes the vector of Euler angles that describe the orientation of {C’} with respect 
LOm pine 

From the definition of € and equations III.9, it can be concluded [Ref. 13] 


that trimming trajectories Po € E correspond to helices: 


0 
a = Lis 
ia : ee 
We 


(2 


—v, cos(¥-) sin(t)ct) 


Wie = | vecos(ye)cos(bet) |, (111.12) 
—v¢ sin(yc) 
where 
w- = desired turn rate, 
D2) —  desited inertial velocity, 
y- = desired flight path angle. 


Thus, the trimming trajectories can be parameterized by the vector 
Ne = [Vey Pe, Fe]? € R?. ier 


In fact, given 7,, we can determine the trimming values for Vo, Qc and Ue which will 
be important later for the controller synthesis. 
Integrating (III.12) we obtain 
- = Ye) os (7).t) 
Pe(t) = =a sin(wt) | - (111.14) 
—v, sin(7<)t 
Equation III.14 indicates that the radius of the helix is 
be = =e (IIT.15) 


and the climb rate is 
h, := —v,sin(+). (111.16) 


Using elementary ideas from differential geometry [Ref. 3], equation III.14 can be 


reparameterized by considering the arclength s defined as 


ssf |ZPold)llat, 


_ / ee (111.17) 
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Although the symbol for arclength is the same as that used in other sections for the 
Laplace transform variable, it is hoped that the context will make the distinction 


clear. Then, using (JII.15),(I1I.16) and (III.17), equation III.14 is written as 


b, cos(b.*) 


* 


Po(s)= | bsin(p-+) | - (I1I.18) 


Two parameters useful in describing the helix, Po(s), are its curvature and its 


torsion. The curvature is defined in terms of the parameters in 7, as 


sar 1).Cos(¢) 


" , (III.19) 
and the torsion as 
here 
as) = ye? 
2 Pee (111.20) 


2. Trimming Trajectory Coordinate System 
Now, let {A} denote a Frenet frame attached to Po(s) [Ref. 3]. The z axis of 
the Frenet frame is termed T’. It is a unit vector tangent to Po(s) at s, and it points 


in the direction of motion on Po(s). Its components are defined as 


Gia 
ms) = “0 
— bts sin(b-) 


berbe cos(Pe=) | > 


Vc 





__ he 
Ve 


— cos(7¥c) sin(y.) 
cos(¥¢) cos(pe) (III.21) 


— sin(y<) 
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The y axis is termed N and is perpendicular to 7’. It is a unit vector pointing toward 


the center of the helix, and it is defined as 


N(sy = S/x(s), 
- cos(v).£) 


sien) | (I11.22) 
0 





I 


The z axis is termed B, and is orthonormal to J’ and N. It is defined according to 
the right hand rule, B(s) = T(s) x N(s): 
sin(7c) sin( ye) 
Bis) = | —sin(y.) cos(bez) | (111.23) 
cos(7Vc) 
The rotation matrix from {A} to {J} can be defined by the vectors, T, N, and B as 


he = (EF NB. (111.24) 


Let 40 be the angular velocity of {A} with respect to {J}, resolved in {A}. It can 
be shown that “2 is equal to [rs 0 «s]’. Using 49 and the Serret-Frenet Theorem 


[Ref. 3], the time rate of change of 4,R is 
WAR) _ a4R), 


dt ds 
D =e) 0 
Sieve 0 —7 | S, 
0 +r 0 
= 4RS(40), (111.25) 


where S(Y) is a skew symmetric matrix defined by [Ref. 10], 


ej Ix. 
0 -—-Ks 0 
a Ks OQ -TS 
0 Ts 0 


Let P(so) denote a point on the trajectory Po € €. Define the error vector 


Mg := 4R(P — P(s0)). (111.26) 


We will call Pe(so) a projection of P onto Po € € when Mg has the following form: 


Mr = {0 y 2)’. 


From the definition of Mg, Pc(so) is also the point on Po nearest to P. (For the 


discussion of when such a point can be uniquely determined, see [Ref. 40]). Let 


Pc (so) be the projection of P onto Po € €. Then, using the definitions of Mg and 
{A}, we obtain 


Therefore, 


aRV 


dt 


7 RRRV, 


nal 
d d, 
io ao q (ak Me), 
g 0 
AR | 0 | +4 RS(4A0)ME+,R| gy |. (III.27) 
0 z 
s 0 


TRAR| 0 | +7 RaRSAOQ)Me +7 RAR | G |, 


Ne 


0 Z 


+ S(40)Mz, 


—YKS 


yTS 


= 2T 1 O y , (II1.28) 


Or 


S l—yk 0 0 
on zr 10 RV. (111.29) 
Za YT O | 


3. Generalized Error Vector 

The design of an integrated guidance and control system for the plant G in 
(111.9) and the set € involves obtaining linear models for G along the trajectories in 
E. These models will necessarily be time-varying in the state space coordinates used 
in (III.9). It turns out, however, that an appropriate coordinate system exists where 
the linearization of the plant G along any trajectory Po € € is time-invariant. This 


coordinate system is discussed next. Let Po, Ac € € be given. Define 


Ve = V—Ve, 
= esos 
r (111.30) 
Pe = | y Z| . 
Ag = oO! [A — Ac], 


which can be interpreted as the generalized error vector between the vehicle state 
and the trajectory in €. Simple physical considerations show that the problem of 
following a trimming trajectory Po € € at a fixed speed Vo is equivalent to driving 
the generalized error vector to zero. 

In order to have a complete description of the error dynamics, we need to derive 
an expression for the time rate of change of the generalized error vector, (III.30). 
Expressions (III.9) and (III.29) suffice to describe the derivatives of the first three 


terms in (III.30), since Ve and Qe are constant along trimming trajectories. The 


T 


fourth term, Ag, is the only term that requires an explanation. Recall, 


Ap =O iw Nele (I11.31) 
Therefore, 
d _ ead d dae 
fhe = O*(Sa- Sr) + £07 [0- dc). (111.32) 
Since 
d 
a () 
dt oe 
and 
d 
ane — OG licr 
we obtain 
d = dager, 
he = 2-Q QoNc + 7 O [A ere (111.33) 


Therefore, the generalized error dynamics are described by the system of equations 


d Roe 
dt Ve = 
d 
a SE 
Ca — 
d = 
dt PE = 
d _ 
dt Ag = 
where 


Fv,(Ve, Qe, Az) ++ Dye Veil at ee ene 
Fon( Ve, Qe, Ag) a5 To, (Ve, In)H( Ve, Nez, U), 


—1 


L—yx 0 0 

ae (111.34) 
ads l 0 ERY, 

ie 0 | 


Oz — NG =O 06Nc +O OAs. 


Fry,(Ve, Qe, Ag) = Fv (Ve + Vo, Qe + Nc, QArz + Ac), 


and similarly for Zy,, Fg, and Zo,. 


An important consequence of this choice of error dynamics is that the lin- 


earization of equations III.34 along any trajectory Po € E is time-invariant. In order 


to derive the linearization of (III.34), we require the following identity. 


78 


Identity 1: Let 1, R be the rotation matriz from { B} to {1}, and let Q satisfy (A) = 
OD. Then for any vector X in R° 


d 

A(BRX) = —BRS(X)O™. (111.35) 
and 

d B ‘ B 7 —1 

(PRX) = S(PRX)Q™. (111.36) 


Proof: First we observe that [Ref. 10]: 


<(bR) =, RS(Q). aes) 


Next, from the definition of a cross-product we get 
Sixes = 2 oy Xx (111.38) 


for any X,Y € R°. Now, using equations III.37, 11.38, and the fact that X is a 


constant vector, we obtain 


d 
dt 


d d 
I ir ee EEL Feed eet 
(RX) = 2GRx +4 REX. 


= pRS(Q)X, 
RSX. (111.39) 


Next, observe that 


er ds... 
d 
= a (BRX)O Q. (111.40) 


Equation III.35 now follows by comparing equations III.39 and III.40. 


To obtain equation III.36 note, 
(7R)* = (BR), 


i? 


thus 


d iB _ _B d B 
“(PR) = —PRO(LR)PR, 
= —§(0)7R. (111.41) 
Therefore, 
oe Bee 
= S(7 RX). (111.42) 


Moreover, using the chain rule, it follows that 


nr! d 
{(rRX) = a(PRX)TA, 


d 
= (7 RX)Q0. (111.43) 


Equation III.36 follows readily from equations III.42 and II1.43. 
Finally, let A = 0, then 


Hea =O ae 


Now, using equations II1I.35 and III.36, we get 


d d 
a (B RX)|a=0 = Pane X )|a=o, 
eG Om 
= S(X). (11.44) 


= 
Since any trajectory Po € € defines a trim condition for equations III.6 and 
III.7, the linearization of (III.6) and (III.7) is naturally time invariant. Introducing 


the following notation, 


AS = SIF) + Ze OHO) 
ao: 
Be = SIKH) 


SO 


the linearization of (III.6) and (III.7) is written as 


“5Ve = A Ve + Ag bE ate A} ONE + By dU, 
£50 = Al5Ve + AM6Ne + APSA g + BobU. (111.45) 


We need only to derive the expressions for £ Pr, and “og to show that they are 


time invariant. 


Let 
=] 
l—ykr 0 0 
aa 2T ne 
YT 0 | 


= | -2z7r(l—ye«)' 1 0]. (111.46) 


Then, from (III.34), 














te = M(Pr)gR(A)V 
=e CY, (111.47) 
Furthermore 
“sp, = — (MARV]|.6Ps + <-[MARVI]|.5Ag + 2o-[MARVI.5V 
dt Ope ae OV a 
(111.48) 
Expanding the first term in (III.48), we obtain 
aM, , 
salMs BRV]|. = op, rc Re, 
ya 9 
— ee BPg, |e 3P a, |e | aRVo, 
00 0 Stew) |, 0 0 0 0 0 
= joe gaer(\my)-* 0 0 aa aes Gunung ey cs 
0 0 0 a-yrltve)" |) gg 0 0 0 
Oy ce 
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00 0 K(ARVc)> 0 0 0 0 0 
= 00 0 0 0.0 =m eye tn 
0 0° 0") = R1 ) Saeano 0 0 0 
O kK 0 
= (2RVe)- | Oman (111.49) 
0 —r 0 


where the subscript c implies that the preceding gradient is being evaluated along the 


trimming trajectory. The subscript + implies the z-component of the vector, and 


T 
Pr = sey lam Ep 
Now, from the definition of Ag, it follows that 
OA 
a Q. (111.50) 
Next, using equation [II.50 and Identity 1, we get 
O O 
——|[MBZRV]|. = ——[gRV))|- 
DAL B l| ae Bre | ’ 
Cs 
= (aa gl ReRV))lc, 
O ON 
= (7R~[pRV]—])|- 
; 5A LB lan) ’ 
= —(7ReRS(V)O"Q)l-, 
= —7RS(Ve), (111.51) 
and 
0 MA A 
By combining equations (III.49 - I11I.52), we obtain 
0 kK OQ 
d 
0 —r Q 
— ARS(Vc)dAzE +4 R6Vg. (111.54) 
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The final term, £éXR, is obtained by observing that, along the trajectories 


Po € €, O71 =0. Using simple algebra, we obtain 


di 
d O d OA 
7) = Oe ——| 6A p. [I1. 
ee . AA di are = a) 


Now, since $Ac = [0 0 J? and from the definition of Q7! (IIL.3) and PR (III.2), 


it follows that 


Ore 
ee C 


I 


Qo, 


d 
7 R= Ac. (111.56) 


Finally, using Identity | and equations II1.50 and II1.56, we obtain: 


6 r\ iaeaed 
(Oa —— ee C) 
aA 2 GSO aig an! PG he) ele 

= S(PRENc)OO|. = S(). (111.57) 


The desired expression for fon now follows from expressions [II.55 and III.57 as 
d 
4 OE = 60 —S(Nc)bAg. (II1.58) 


Summarizing, the linearization of equations [11.34 along any trajectory Po € € 


is time-invariant and given by 


£6V_ = AV6Ve + AY6Qe + AV 5Ag + BySU, 
260, = AS6Vp + AX5Ng + AP6AR + Bo SU, 


0 Kc OQ 
CG in OP ide ve). | 0) 0 —7 | 6Pz, (111.59) 
0 -r 0 
— ARS(Vo)éAg +4 ROVz, 


| “éAr = 60g —S(Nc)éAg, 


where (ARVc), is the z-component of the vector ARVco. In the next section, the 


symbol G, will denote the set of all linear plants G, (n,.) associated with the set €. 
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D. TRAJECTORY TRACKING CONTROL SYSTEM DE- 
SIGN AND IMPLEMENTATION 
1. Linear Controller Design 


In the previous section, we have shown that the linearization of the nonlin- 
ear system Gg about any trajectory in € results in a time-invariant plant G (7). 
Therefore, associated with the set € there is a family of linear plants, G, . These can 
now be used to synthesize a tracking (possibly gain-scheduled) controller C , which is 
designed to operate over all the trajectories in €. 

A common approach to the development of such a controller C requires de- 
signing a family of linear controllers for a finite number of linear plants in G; , and 
then interpolating between these controllers to achieve adequate performance for all 
the linearized plants in G; . During real time operation, the controller parameters are 
updated as functions of a gain-scheduling variable q = h(V,2, A, P,U,n-.), where h is a 
C’ function. For example, a typical gain-scheduling variable, for the case of aircraft, 
is dynamic pressure. In that case, q = $p||V||?, where p represents air density and is 
itself a function of aircraft height. 

In what follows, we restrict ourselves to the idealized case where the description 
of each controller for each plant G; (n.) is available [Ref. 39]. Therefore, we assume 


that the first design step produces the set 


C) = {C; (gainqe aa h(Vc, Qo, Ac, Po, Uh ae 


given by (see Figure 33): 


6E = [dv dy Sz]? —[dv, dy, 52,]*, 
£5X eI — Aa(gce)dXe1 ate Ba(q-)[6V7 508 Oa 
+B.2(q-)6X c2 ae B.3(q-)6E, 


Ci (qe) = (111.60) 


d ie 
70% 2 = 6B 


? 


a Cori(qe)bXe1 a> Da(qe)[6V7 én? bAR]* 
+D.2(9c)6Xc2 + De3(qe JOE, 
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where 6X., € R™, 6X.2 € R™ and m = dim(U). The vector 6X2 represents the 
integrator states of the controller C; (q-). The vector 6X,, represents the remaining 
states of the controller C; (q-). The velocity and positon error is 


Ht 7 





by my 9 OF] 6V 0 0 0 
by = 0 0 0 SO | + |e0eel. GO } OPE, 
52 0 0 01] | éAg Dai 
6V 
=: C,| 62 | +C.6Pz, (111.61) 
6\g 


where év,, dy, and éz, are introduced to determine how fast the error states, dv, dy, 


and 6z, go to zero. We further assume that the parameters of the controller are C! 


functions of qe. 







OVE SO oe 





[duc byc bZ¢ 





[dv dy dz]? 


_ Aei(Ge) Ba(qe) Be2(qc) Be3(qe) 
Ka) = | gates Da(qe) Dea(ae) patie 


Figure 33. Set of Linear Controllers 
The structure of the controller C; (q,) has the following important feature. 


Suppose the closed-loop system consisting of G (q-) and C; (q-) given by equations 
11.59 and III.60 is asymptotically stable. Then, for a given q,, the controller C; (q,) 
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will ensure zero steady-state error to a step input for the variables in 6F. This includes 
errors in the vehicle’s inertial velocity v, and in the deviations from Pc, z and y. Zero 
steady state errors are achieved by integrating 6. This structure is typical of tracking 
controllers, since they are designed to drive errors between step changes in reference 
commands and the corresponding plant outputs to zero in steady state. Notice that 


the block K(q_) (see Figure 33) may itself contain additional integrators. 


2. Gain Scheduled Controller Design 

Next, the family of linear controllers C; (q.) must be implemented on the non- 
linear plant G defined in Section B. This problem has been addressed in [Ref. 27] 
for the general class of nonlinear plants and for tracking controllers with the same 
structure as C; (q-). In [Ref. 27], the authors formulated a so-called controller imple- 
mentation problem which will be repeated here for the problem at hand. 

Let T(G: (qc), Ci (d-)) be the closed-loop linear system that results from con- 
necting C; (q-) to G (q-), and denote by T(G (q-), Ci (q-)) the corresponding matrix 
transfer function. Let T7(G , C )(q-) be the nonlinear closed-loop system that consists 
of C and G, and let 7)(G , C )(q.) denote its linearization about Po € €, and denote 
by Ti(G , C )(q-) the corresponding matrix transfer function. With this notation, the 
controller implementation problem applied to the integrated guidance and control 


problem considered in this chapter can be stated as follows: 


Controller implementation problem: “Find a gain scheduled controller 
C such that for each trajectory Po € € 


1. the feedback systems 7)(G , C )(q-) and T(G (qc), Ci (q-)) have the same 


closed-loop eigenvalues. 
2. The closed-loop transfer functions T;(G , C )(q,) and T(G: (qc), Ci (qc)) are 


"equal.” 


Given the set C; of linear controllers for the set G, of linearized plant models, 
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we propose the following structure for the gain scheduled controller C (see Figure 34): 


[Oy 2] = 7R(P— Pe(so)), 


— a: 
4Xa = Aa(qg)Xat+Bal(g[Zv? £07 (Q-O- Ag)? I]? 
C i= +B.2(q)E 4 B.3(q). 8, (111.62) 


Fe 
£X2 = Cal(q)Xe + Da(q)[EV7 
+D.2(q)E + D.s(q) SE, 
U = X 2: 


LOT ( = Om Neale 


1 


Recall, Pe(so) is the projection of P onto the helix Pe € €. Comparison of Figure 33 
and Figure 34 indicates that the structure of the gain scheduled controller is easily 


obtained from that of the linear controllers. 





Figure 34. Gain Scheduled Controller 


We now make the following assumptions: 
Au, (Dyin Coy) alien) —=“obbaa aa) 
A2. The matrix 
sl -Aa(q) Be2(q) 
| ni 4)) D.2(q) 





S71 


has full rank at s = 0 for each Po € E. 
A3. The matrix pair (Aq,C.1) is observable. 


Assumption Al implies that the number of integrators is equal to the number of 
control inputs. This is necessary if the controller is to provide independent control of 
the errors & using the control inputs U. Assumption A2 implies that the realization 
(Aq, B.2,Ca,DP-2) has no transmission zeroes at the origin. Finally, assumption A3 
guarantees that the state X,) is zero along the trajectories in €. 


The main result of this section is stated next. 


Theorem D.1 Suppose assumptions Al, A2 hold. Then the gain scheduled controller 
C given by equations III.62 solves the controller implementation problem, 1.e., for each 
Po € € the following properties hold: 


1. The feedback systems 7)(G , C )(q-) and T(G (q-), Ci (q-)) have the same 


closed-loop eigenvalues. 


2. The closed-loop matrix transfer functions 7;(G , C )(q-)(s) and T(G (qc), Cr (de))(s) 


are equal. 


Proof: In the proof, we set the controller matrices Du, D.3 to zero. This does not 
change the results but considerably simplifies the algebra. 
Let Po € € be given. Let 


Ae Ae AN 
A = A} AC Ae ) 
0 IT —§(M) 
A, = | AR 0 -ARS(Vo), |. 
O K 0 
A3 = iia 0 0 car ’ 
0 —r Q 
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By 


b — Bo 5 
0 
(ARVc)z := x-component of ARV 
bVe 
One n= 6Qr ’ 
bArR 


Then, (111.59) is expressed as 


£dY = Aj,6T + BSU, 
Gi (Ne) = (111.63) 
£§P = A2dT + AP. 
The error signal is expressed as 
0f = C61 +C,6P, (111.64) 


which, when substituted into the linear controller (III.60), results in 


en = Aal(qe)dX ci =f B26 X 2 a [Bei (qe) az Be3(qe)Ci| éY 


@(q.) = et (111.65) 
£hX.2 — C,67 + C26 P, 
1 Cor (de )EX 1 Deo (qe)6Xc2- 


Consider the feedback interconnection of the linear plant (III.63) and linear 


controller (III.65). The state matrix F' of this feedback system has the following form: 


A BD 3Ci BCa BD 2 


Ay A3 0 0 
F:= | (111.66) 
Ba Be3Cy Be3C2 Aa Be2 
Cy C2 0 0 


Next, we linearize the feedback interconnection of the plant G and the con- 


troller C . However, in order to do that, first we must determine the values of the 
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controller states X., and X,2 along the trajectory Po € €. From equation II1.62, we 


obtain: 
d - d T d Te Oa!) Tag 
Gxee = Acrlqe)Xac t+ BalgellZVe GM (Ne — Qe'Ac)") 
d 
+B .2(q.) EL al Bela) E, 
d : 
Gyo ede = Cal@oeer + Deo(qe)£, 
Ui. = Ae (111.67) 
Notice, since along Po € E: 
eae a) 
Vo = “Colstanse, 
eee re olisvantr 
Ng -OQg'Ac = 0, 
X 5. a Ua constants (111.68) 


we get 


ae a c NG ) 
dt le A 1 le 


0 c= Creer 
Now, using Assumption A3, we conclude that 
SG — 0. (II1.69) 


In order to compute the linearization of the feedback interconnection of G 
and C (7T(G,C )) along Po € E, first observe that F(G ,C ) is equal to the feedback 
interconnection of Gg and the system consisting of A(q) and the integrator X.2. The 
linearization of Gg along Po € E is given by (III.59). The linearization of the system 


consisting of A(q) and X.2 can be obtained using the trimming values of X., and X-2. 
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Define 














d . _ 
gl = [SVE 498 (@-o-%ey" | 
T 
ao := | x3 xe <0 ae iE de 
then, along Po 
i 
Tole Lo UT 0 0 0 deo | (111.70) 
and the linearization of (III.67) can be written as 
d., _ HAa(ge)Xer.), _, ABalae)£0, , O(Bea(ae)B) 
ee as CU TS 
O(Bes( qc) E) 
7a — |a03 
Ow 
a. _ ACalge)Xer), _, (Deal 4)E) 
a SC 
_ O(X-2,) 
a A on Tere) 
Consider the expansion of the first term in (III.71), 
O(Aei(qc) Xr.) as O(Aa(qe) Xr.) Of Acr(de) Ker) | 4 QA lac) Xare) 
Ow “3 oe Oe ° (£6) 
O(Aer(Ge)Xe1.) O( Aci (Ge) Xt.) O(Aei(qe)Xer.) 
+ OE | + aE Wo ala Aq. Ge 
= Ari(ge Xe. (111.72) 


Since X,1, 1s zero along the trajectory, their are no extra terms due to the scheduling 
parameter, g-. Similar results are obtained for the terms B.4(q-), Ca(qc), and De2(qc). 


Noting that 
d if 


pees 
bo = | GOVE FON FONE | > 
the linearization of (III.67) has the following form: 
Cae ee ee 2. d 
ce = AadXea. + Bala oVe qe apoE! + B2dk + Baa ok, 
d 
Gy oke2 = C16X1. + D2dbdL, 
OU = OLE (111.73) 


9] 


The state matrix M of 7)(G , C ) is calculated as 




















A, 0 0 B 
A, Az 0 0 
Mi i A, 0 ia B 
[Be B.3| + Beo\Cy Co] Ac (Ba B.3| 
C, Cy A, Az Ci Cy 0 
De2 [Cy Co] Cr 0 


The proof of the first part of the theorem can now be shown by following the proof 
of Theorem 4.1 in [Ref. 27]. From Assumption A2, it follows that the matrix 


Ac B.2 
Cy Deo 








is invertible. Therefore, we will use this to define a very useful similarity transforma- 


tion to show that F’ and M have the same the eigenvalues. Let 


—1 














Aa Be er ale 
= ; (III.74) 
Cai De Z W 
and set 
I eat) 
I O 
ee Ba Bes x ie 
= pe (IIL 
| I 0 
Bo | Ba Be ZW 
| Ci; Co 
Expression II].74 implies 
Aa + B.2 — ie 
Ce Yo Deo a de 
Ag xg ae Be = 0, 
CaXx + D2 = 0. (111.76) 


o2 


A simple computation using (III.76) and (III.75) verifies that 





, 0 0 
1 0 
ica |B. Ba | x y |, (111.77) 
CC 
0 ZW 


and that F =7TMT—'. Thus, F and M have the same the eigenvalues. 
In order to prove the second part of the theorem, it will suffice to show that 


the linear controllers (III.73) and (III.60) have the same transfer function from 
T T 
—— | by dy 5: | —| bv. bY, $2. | : 
and 
T 
— 6bVi §6QF SAT (111.78) 
to 6U. A direct calculation of the Laplace transforms show that 


U(s) = Coi(sl — Aye {Ba 6V7(s) 607 (s) 6A7(s) | 


Bok T 
a 5 Bab(s)| ae 6V7(s) 5N7(s) SAT(s) 
_ 
— + DsE(s), (IIT.79) 
for both controllers. a 


Thus, the eigenvalues of the linearizations along each trajectory in € are pre- 
served. Furthermore, the input-output behavior of the linearized operators is pre- 
served in a well-defined sense. The reader is referred to [Ref. 27] for a complete 


discussion on approximations to this method that avoid using pure differentiation. 


3. Implementation Procedure 
The Theorem D.1 can be used as follows: first, determine the dynamics of 


the vehicle G and the set of trimming trajectories € the vehicle is required to track. 
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This set is parameterized by the range of trimming velocities (v,), desired flight 
path angles (y.), and desired heading rates (~),). Recall, these variable constitute 
a gain-scheduling vector gq. = [Vc Yc pe]. Next, rewrite the vehicle dynamics 
using generalized error coordinates Vg, Qe, Ag, and Pr to obtain the vehicle’s error 
dynamics Gg. Linearize Gg about a finite number of trajectories in €. Use these linear 
models to design a finite number (k) of trajectory tracking controllers {C; ;,2 = 1, k}. 
Gain schedule {C; ;,2 = 1,k} utilizing a favorite interpolation or gain-scheduling 
technique to obtain the linear gain-scheduled controller C; (q.). Now, implement this 
controller on the nonlinear plant G according to the expression III.62. 


It is worth emphasizing the following important properties of the controller C : 


e The result in Theorem D.1 holds for all trajectories in €. 


e The structure of the controller C is easily obtained from that of the linear 
controllers. 


e Since all the closed-loop transfer functions of the local linearizations are pre- 
served, at the level of local linear analysis, the controller does not introduce 
any additional noise amplification despite the presence of a differentiation op- 
erator. 


e Along trajectories Po € €, X.2, = U, and X.., = 0. Therefore, the trimming 
values of the control inputs are naturally provided by the integrator block with 
state X.2,. 


e The integrators X,2 are directly at the input of the plant, which makes it 
straightforward to implement anti-windup schemes. This becomes necessary 
in applications where the input U is hard limited due to actuator saturation, 
for example. 


e ‘The inputs to the controller, Po and ee can be computed directly from the 
vector qe. 


e The trim values Vo,Qc, and U, are not required in the controller implemen- 
tation. 


e Along the trajectories Po € €, the controller guarantees that the steady state 
value of error vector F is zero, which follows from the fact that the controller 
solves the controller implementation problem. This is in sharp contrast to 
standard LOS guidance schemes. 
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E. EXAMPLE 
In this section we apply the methodology developed in Section D to the design 


and implementation of an integrated guidance and control system for a fixed wing 
unmanned air vehicle. The vehicle’s stability and control derivatives are the subject 


of a section of chapter IV. Using the notation in Section B, the Frog dynamics have 


the following form: 


$V =Fy(V,0,A) + Ty(V,Q)H(V,Q, U), 
ee Va ee fal V, 9) VQ, U). 

ee ; a( ) + ZolV, OVH( (111.80) 
ae eV. 
eA =O? 


Following the development in Section B, the set of trimming trajectories € for 
the vehicle is defined as follows: 
ie 
Ac 


+ Pe = IVE. 

£ ho = Qc Nc, 

Fy (Vo, Qc,Ac) + Zv(Ve, Qe )JH(Ve, Xe, U.) = 0, 
0 


POWVeianiclet i olVa, ic \i{Ve,2c,U0.) = U, 
(111.81) 








= 


where Pc and Ac can be computed using the scheduling vector gq = |v. 7% pe]. 
Now, given [v. 7 be]? and 8, = 0, we can solve for Vo,Nc, U., and Ac: 
Fvy(Voe,Qe, Ac) + Iv(Ve,2c)H(Ve, Ic, U-) = 0, 
Fo(Ve, Ne, Ac) + Za(Ve,Qc)H(Ve, Nc, U.) = 0, 
O51N¢ —Ac =0 
Vol] — ve = 0, 
B.— sin“! = = 0, 


t 


Ye — [0 1 0] arg(y R7 R) = 0, (111.82) 


where the arg function extracts the angles X from the rotation matrix R(X): X = 


arg(R(X)). 
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Using the solution to equations III.82, the linear model for the vehicle rep- 
resented by equations III.59 was obtained along a trajectory characterized by the 
velocity of 73 feet per second, flight path angle of zero, and heading rate of 10 degrees 
per second. This model was used to design a linear trajectory tracking controller for 


the vehicle. 


1. Design Requirements and Linear Controller De- 
sign 


Design requirements for the linear trajectory tracking controller included 


e Zero Steady State Error: Achieve zero steady state tracking errors of all 
trajectories in €. Achieve zero steady state tracking error of indicated airspeed 
while on any trajectory in €. 


e Bandwidth Requirements: The command-loop bandwidth for each com- 
mand channel should be no greater than | radian per second and no less than 
1/10 radian per second; the control-loop bandwidth should not exceed 12 ra- 
dians per second for the elevator, aileron and rudder loops, and 5 radians per 
second for the throttle loop. These numbers represent 50% of the correspond- 
ing actuator bandwidths, and shall ensure that the actuators are not driven 
beyond their linear operating range. 


e Closed Loop Damping and Stability Margins: The dominant closed- 
loop eigenvalues should have a damping ratio of at least 0.5. Simultaneous 
gain and phase margins of 6db and 45 degrees in each control loop must be 
achieved. 


The methodology selected for linear control system design was H. synthe- 
sis [Ref. 12]. This method rests on a firm theoretical basis, and leads naturally to 
an interpretation of control design specifications in the frequency domain. Further- 
more, it provides clear guidelines for the design of controllers so as to achieve robust 
performance in the presence of plant uncertainty. The basic steps in the controller- 
design procedure, including the development of the synthesis model, were done us- 
ing the approach described in [Ref. 24]. This approach provides an intuitive and 


straightforward way for converting the design requirements into the weights for the 
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H.. synthesis model. Consider Figure 35. Here C; is the controller to be designed, 


G, is the linear model of the vehicle. 





Figure 35. Synthesis Model 


In Figure 35, the vector of exogenous inputs w represents the commanded 
inputs. The vector y, represents lateral and vertical displacement states of the linear 
model as well as the vehicle’s velocity. The regulated output z includes the outputs 


of the weighting matrices W; and W2. These matrices had the following form: 


0.0 
W, = Oe Or); 
ce 
cantly © 
W, = 0 <a (111.83) 
0 0 & 


where the constants c;, 2 = 1,6 were used as the design knobs adjusted to meet the 


closed-loop tracking, damping, control, and command loop bandwidth requirements. 


om 


Values of W, and W 2 were iterated on, and used to obtain a linear trajectory 


tracking controller: 


é6E = [dv dy bz]? —[du, by bz,]7, 
£OX 1 = AadXe a B16 Xe; 
LX o2 = Of, 
— 6U = ®{Ca6Xa + Da[6V% 627 AR)? + Dood Xc2 + DexdE}, 


Ci (q, ) = 


where the H,, state-feedback gain is K = [C. Da De Des]. The feedback system 
consisting of the plant G, and the controller C; was found to meet all the design 
specifications given earlier in this section. Since the control surface effectiveness is 
proportional to the dynamic pressure, the controller C; was gain-scheduled on @. 


The value go represents the nominal value of g. 


2. Implementation and Simulation Results 
Using the formulae provided in Section D, the family of linear gain-scheduled 


controllers C; (¢) was implemented on the nonlinear plant G as follows: 


(Oy 2) 7 1 Foca) 


E se lv — Usd Al 
Ci=4 Xa =A Aq + Bak, 
re =a {Ca Xe =p DalZV7 £QT ) ~ O-(A)ZAc]? Ss DatE oe DoE}, 


U =" 5) 
We emphasize that the implementation equations for the controller C do not require 
the computation of Qc, U. and Vo. Moreover, since Nc = [00 belt, the controller 
must only be provided with wv, and Po when steering the aircraft along the trajectory. 
These are the critical advantages of the proposed methodology. The acceleration 
term £V can be computed using onboard sensors without resorting to differentiation. 
Therefore, the only term which could not be computed directly was £0. In this case, 


the differentiation operator = was replaced by a causal operator with the transfer 


function — 5 [Ref. 27]. 
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The trajectory tracking controller developed above was tested using a number 
of trajectories. One such trajectory consisted of a straight line transitioning into a 
helix shown in Figure 36. This trajectory is characterized by a cruise velocity of 73 
feet per second. Initially, the trajectory is aligned with the inertial x-axis. After 
proceeding along the x axis for 3000 feet, the trajectory turns into a helix with a 
radius of 1000 feet and climb angle of 5 degrees. Consider Figure 37, which shows 
the time history of the position error, bank and pitch angles, and indicated airspeed 
along the trajectory. Clearly, the controller drives the vehicle along this trajectory 


with zero steady state position errors while maintaining 73 fps indicated airspeed. 
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Figure 36. The trajectory tracked in simulation. 


F. CONCLUSIONS 

A new method was introduced for designing and implementing integrated guid- 
ance and control systems for autonomous vehicles. The starting point is a family of 
linear controllers with integral action designed for linearizations of the nonlinear equa- 
tions of motion described in an appropriate state space. Based on this family, the 


method produces a gain scheduled controller that preserves the input-output proper- 
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Figure 37. Time history of position errors, Euler angles and airspeed along the 
trajectory. 


ties of the original linear closed-loop systems as well as the closed-loop eigenvalues. 
The key feature of the method is the ability to automatically reconfigure the control 
inputs of the vehicle to provide for proper control action as the body tracks an inertial 
trajectory in free space while maintaining constant airspeed. The method is simple 
to apply and leads to a nonlinear controller with a structure similar to that of the 


original linear design. 
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IV. RAPID PROTOTYPING SYSTEM FOR 
FLIGHT TEST OF AN UAV 


A. INTRODUCTION 
This chapter describes the development of a rapid prototyping system for flight 


testing of guidance, navigation, and control algorithms for unmanned air vehicles. The 
system affords a small team the ability to take a new concept in guidance, navigation, 
and control from initial conception to flight test. In order to do this, a number 
of engineering problems had to be overcome addressing a gamut of issues including 
weight, power, portability, risk, electronic interference, vibration, manpower, etc. The 
main contribution of the project is the proof of concept flight test demonstration of a 
new integrated guidance and control algorithm (see Chapter II). The success of this 
endeavor had a synergistic effect on several other projects, paving the way for joint 
ventures in voice controlled flight, coastal mapping, and autonomous landing. The 
project is viewed as the foundation of a long-term investment in innovative unmanned 
air vehicle applications leveraging, in part, the operational experience of the officer 
students in the avionics curriculum. 

Testing of a new algorithm, sensor package, vehicle, etc., requires expertise 
from many branches of the engineering sciences, especially aeronautic, electrical, and 
computer. It is potentially costly and time consuming, as well as having the potential 
for catastrophic failure. When successfully done, however, it provides developmental 
information, insight and data that are unavailable from other sources. All of our 
theoretical and numerical results must be verified by some form of experiment, and 
flight testing is often the best way to do this. 

The chapter begins with a conceptual discussion of the Rapid Flight Test Pro- 
totyping System (RFTPS). Motivation for its development is addressed. Next, a 
description of the hardware components is given. The main contribution of the chap- 


ter is discussed in the last section, where an application of the RF TPS to the problem 
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of integrated guidance and control introduced in Chapter III is presented. The full 
capabilities of the RFTPS are demonstrated when this novel guidance algorithm is 


taken from theoretical development to flight test. 


B. SYSTEM DESCRIPTION 


The RFTPS consists of a test bed unmanned air vehicle equipped with a com- 
plete avionics suite necessary for an autonomous flight as shown in Figure 38, and of 
a ground station responsible for flight control of the UAV and flight data collection as 
shown in Figure 39. A functional block diagram of the RF TPS is shown in Figure 40. 
The key goal was to use off-the shelf-technology as much as possible, thus exploiting 
the economy of scale of a number of commercial industries. Furthermore, if the UAV 
development program is to span many years, and to draw on the talents of the officer 
students in the future, the RF TPS had to emphasize high level algorithm design. Low 
level code and device driver generation is kept to a minimum with the vast majority 
of the code ’writing” being done via autocode tools. The system architecture is open, 
providing the ability to add, remove or change real time input/output (I/O). Compu- 
tational power can be increased as mission requirements dictate. The telemetry links 
are secure, yet low power and unobtrusive to the public, not requiring advance per- 
mission for use, special frequencies from a government authority, or special airspace. 
The onboard components are light weight and low power, allowing for the inclusion 


of additional payload. 


Pe RFTPS Capabilities 
The RFTPS developed provides the following capabilities. 


e Within the RFTPS environment, one can synthesize, analyze and simulate 
guidance, navigation, control, and mission management algorithms using a 
high level development language. The same code that ran the simulation, flies 
the vehicle. 


e Algorithms are seamlessly moved from the high level design and simulation 
environment to the real time processor. 
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Figure 38. Frog: The unmanned air vehicle at the Naval Postgraduate School. 


The RFTPS utilizes industry standard I/O including digital to analog, analog 
to digital, serial, and pulse width modulation capabilities. 


The RFTP5 is portable, easily fitting in a car. In general, testing will occur 
at fields away from the immediate vicinity of the Naval Postgraduate School. 


The unmanned air vehicle can be flown manually, autonomously, or using a 
combination of the two. For instance, automatic control of the lateral axis can 
be tested while the elevator and throttle are controlled manually. 


All I/O and internal algorithm variables can be monitored, collected and an- 
alyzed within the RFTPS environment. 


2. Cost, Safety and Other Considerations 


Cost and risk are two leading, and at times competing, concerns that had to 


be effectively handled. Since initial testing is to occur within line of sight at all times, 


a pulse width modulated (PWM) remote control system manufactured by Futaba was 


chosen. Testing of a new control algorithm is similar to handing over control of the 
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Figure 39. The base station of the RFTPS in use at the airfield. 


aircraft to a student pilot. The algorithm should have full freedom to perform, yet 
adequate safeguards must exist 1n case it fails. With some modifications, the extensive 
master-slave flight training capabilities built in to the existing RC transmitters were 
exploited. A significant portion of the cost of the RFTPS resides in the real time 
processor, I/O board and modules, and in the host computer. Additionally, while 
compact, the weight and power requirements of these components are significant 
when compared to onboard power and payload available. In order to gain additional 
payload, and in order to manage the risk associated with the loss of an expensive 
computer package, the real time controller was kept on the ground. Sensor and 
control links to the real time controller were bridged via RF components described 


later. 
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Figure 40. RFTPS Hardware Architecture 


3. Components 

The centerpiece of the RFTPS ground station is the AC100/C30 system from 
Integrated Systems Incorporated. The key feature of this product is its autocode 
tools. With a relatively short time available for research by the officer students, em- 
phasis had to be shifted from code writing, debugging and maintenance to algorithm 
development. AC100/C30 utilizes ”Xmath/SystemBuild”, a graphical programming 
environment that uses a high level block diagram paradigm for modeling of linear 
and nonlinear systems. Within the ”Xmath/SystemBuild” environment, the algo- 


rithm can be built, simulated, tested, and debugged. Real-time code can then be 
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generated for execution on the real time processor. 

Currently, the ”Xmath/SystemBuild” environment resides on a Sun worksta- 
tion. PC’s running Windows NT are also supported. This allows for a significant 
downsizing of the ground station should resources become available in the future. 
Communication with the real time processor is via an Ethernet bus using TCP/IP 
protocol. AC100/C30 provides excellent animation tools for building graphical user 
interfaces (GUI). Through the appropriate design of these interfaces, the flight test 
team can monitor, modify ,and control the actions of the real time processor. The 
GUI resides on the workstation. Communication between the workstation and the 
real time processor is via the Ethernet connection, and is managed by a host PC. Ad- 
ditionally, the host PC provides power to the real time processor, as well as providing 
utilities for compiling, linking, and downloading the C-code. 

The I/O consists of four multi-mode, bi-directional, serial ports utilizing RS- 
232 protocol, a 16 channel pulse width modulation port capable of measuring up to 
sixteen PWM signals or generating up to six PWM signals, and a six channel digital- 
to-analog converter. The I/O modules are hosted by the same PC that holds the 
real time processor. The real time processor is a single Texas Instruments Digital 
Signal Processor (TM5320C30). The capability exists for upgrading the processor, 
or running multiple processors, to meet computational demands of future projects. 

The control configuration of the air vehicle is conventional with three indepen- 
dent surfaces (elevator, aileron, rudder) and a throttle. Manual control is provided 
via a Futaba, dual conversion, PWM transmitter utilizing the portion of the radio 
spectrum reserved for Radio Controlled (RC) flight, 72.030 MHz to 72.990 MHz. 
Precautions entail a search of the electronic spectrum utilizing a hand held spectrum 
analyzer, as well as standard procedures employed by RC hobbyists to avoid two in- 
dividuals selecting the same frequency locally. Built in capabilities of the transmitter 
include the ability to transmit one or more signals from a slave transmitter. The 


slave transmitter 1s a modified Futaba transmitter where the manual control effectors 
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have been replaced by a direct connection to the digital-to-analog I/O module. In 
this way, an exogenous source (RF TPS) can be given control of one, some, or all of 
the control actuators of the aircraft using the same RC link currently controlling the 
aircraft. 

The sensor suite onboard the air vehicle consists of an inertial measurement 
unit (IMU) with a three axis rate gyro, three axis accelerometer, magnetic heading 
indicator, two axis pendulum, phase differential GPS receiver, elevator, aileron, and 
rudder actuator position sensors, and angle of attack, side-slip angle, pitot-static, and 
static pressure air data sensors. A four channel analog-to-digital converter is used to 
capture any four of the following sensors: elevator, aileron, rudder actuator position, 
angle of attack, side-slip angle, dynamic pressure, or static pressure. 

Communication between the sensors and the real time processor is via a serial 
link. Low power, matched, spread spectrum RF links provide up to 115 Kbaud rates 
at over 10 miles range. They require no license, and can be used anywhere in the 
United States. Additionally, the onboard GPS unit maintains contact via a serial link 
with a GPS receiver on the ground, which provides differential corrections. 

Power for the onboard avionics is supplied from a lithium-ion battery. The 
battery provides 6 hours of continuous use. The power budget of the onboard com- 


ponents is shown in ‘Table III. 


[Voltage [__Current_| Power | 
[_IMU__| 24 Volts [ 400 Milliamps | 9.6 Milliwatts | 
[__GPS__[12 Volts [ 250 Milliamps | 3.0 Milliwatts_| 
[ Telemetry 1 | 24 Volts [ 300 Milliamps | 7.2 Milliwatts_| 
| iwatis_| 
| | 

| 





Telemetry 2 | 24 Volts | 300 Milliamps | 7.2 Milliwatts 
Talk | )~|~”*~*~*~*«*~SC( Miilliwatts 
[Battery [S200 mW 


Table III. Power Budget of Onboard Components 
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C. TRAJECTORY CONTROL: AN APPLICATION 

The intent of this application was to demonstrate the utility of the RFTPS by 
flight testing a new integrated guidance and control algorithm that was shown to have 
good performance and robustness properties both theoretically and in simulation (see 
II). This project was chosen because, on one hand, it is a totally unique application 
in terms of the guidance and control algorithms implemented. On the other hand, 
to implement this algorithm, most of the steps required of any application involving 
autonomous flight of an air vehicle must be accomplished. If done correctly, the 
foundation will be laid for follow on projects and joint ventures utilizing unmanned 
air vehicles. 

Generic tasks fundamental to guidance, navigation and control algorithm de- 
velopment were accomplished first. To begin with, a high fidelity model of the test 
bed vehicle, nicknamed Frog, was developed. This involved a complete lateral and 
longitudinal parameter identification of Frog’s stability and control derivatives. This 
lead to the development of a six degree of freedom, nonlinear simulation. In order 
to manage the risk associated with the autonomous flight, it was decided to add an 
onboard inner-loop controller. The inner-loop controller modifies the vehicle dynam- 
ics such that displacements of the primary longitudinal and lateral control effectors 
correspond to various trimming trajectories, as defined in Chapter IIJ., The inner-loop 
controller selected was a commercial autopilot whose control laws had to be identified 
as well. Additionally, accurate and timely calibration of the I/O needed to become 
an integral part of the RFTPS to account for changing environmental conditions, 
battery levels, etc. 

This section begins with a brief background of the parameter estimation algo- 
rithm used to identify a model of Frog. Next, the experimental testing completed to 
obtain data used in the parameter identification process is summarized. Flight test 
data is compared with the simulation results to demonstrate the efficacy of the ef- 


fort. Then, a linear controller design based on modeling results obtained is discussed. 
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Next, a nonlinear controller implemention based on the theory developed in Chapter 
III is presented. Finally, flight tests conducted using the RFTPS to track an inter- 
esting representative trajectory are discussed. Data from flight tests are presented 
and contrasted with results from simulations. Along the way, implementation issues 


germane to the flight testing process are explained. 


1. Theoretical Background for Model Identification 

Parameter identification of the test bed vehicle and inner-loop autopilot was 
accomplished based on a series of test flights and lab experiments. The Maximum 
Likelihood Method was chosen to identify the parameters of the model. The back- 
ground is extracted primarily from [Ref. 18]. Assume that the model to be identified 


has the following form 


z = A(p)z + Bip)u, 
y = C(pjz+ D(p)u, 


oe stale Vector, 
u = input vector, 
y = otput vector, 
Di DakaimMeter VECLor. 


The problem at hand is to estimate the parameter vector, p, given a data set of 
measured inputs, u, and outputs, y. To that end, denote the estimated output based 


on the parameter vector p as j(p), and use it to define a cost function 


§(p) — yll3- 








J(p) = 


The method seeks to minimize the cost, J, by varying the parameter vector p. Given 
an initial guess of the parameter vector po, the Jacobean of the cost function J is 


calculated by numerically perturbing the parameter vector po 


OJ _—-y (Po + Op) — y(Po) 


Op ép 
=: Ho. 
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Then, a linear approximation for y(p + ép) is 
y(p + 6p) = y(po) + Hoédp. 


The cost J can now be expressed as a function of the perturbation of the parameter 


vector as 


o(p + 6p) — yo, (IV.1) 
(Po) + Hodp — yl, 
\3(po) — yll3 + 29(po)* Hodp — 26p" Hey + dp" Hodp. 


J (dp) 


In order to determine a 6p which minimizes equation IV.1, set 


dJ : 
Cn 2H6 9(po) — 2Hg y + 2Hg Hoép, 
= 0, 
from which it is evident that 
Hy Hodp = Hg (y — (po) (IV.2) 


must hold. 

A common problem encountered in solving a parameter identification problem 
is the presence of redundant parameters in the model. This results in the matrix 
H? Ho being singular, and presents a problem if a solution to equation IV.2 is sought 
in terms of ép. In order to avoid this singularity, the following problem is solved 


instead, 
(Hj Ho + ul)dép = Hi (y —G(po)). (IV.3) 


With the descent direction dp known, a line search is used to calculate the step size. 
This method is guaranteed to converge to a minimum for the cost function J, although 
not necessarily a global minimum [Ref. 18]. 

The appealing feature of this method for the problem at hand is that it can 


exploit a great deal information known about the model. For instance, the order of the 
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plant and most of the structure in the system matrix is known, as is a reasonable range 


of values for the parameter vector p based on empirical methods (USAF DATCOM). 


D Vehicle Model Identification 


The conventional configuration of Frog suggested that the parameter identifi- 
cation problem could be decoupled into longitudinal and lateral dynamics. Maximum 
likelihood parameter identification was used to refine existing analytic estimates of 
the longitudinal stability and control derivatives [Ref. 34]. The results for a few key 
longitudinal stability derivatives are compared to analytic estimates in Table IV. The 
primary difference was in the estimate of the elevator effectiveness, which was less 


than half as effective has previously thought. 


[Derivative [Analytic | Experimental | %2 | 
[ o, | 490 | 409 [513] 
[Cm | -417_| 557 | -25.1 | 
[C. | ti [301 —*| 186 
[ Cm, | 1.62 [1.05 | 543 | 


Table IV. Comparison of selected longitudinal derivatives. 










Figure 41 shows the response, in simulation, of a longitudinal model of Frog 
using these derivatives to an elevator doublet. Comparisons are made to the results 
from flight test. 

Next, the process was repeated for the lateral axis. Aileron and rudder dou- 
blets were executed while aileron and rudder position, roll and yaw rates, and side-slip 
angle were measured. Results, for a few key lateral stability derivatives, are compared 
to analytic estimates in Table V. The biggest difference was the value of Cy,, which 
increased by over 300 percent. This doubled the estimate of the natural frequency of 
the dutch roll mode. Figure 42 shows the response, in simulation, of a lateral model 
of Frog using these derivatives to an aileron doublet. Comparison to the results from 


flight test are also shown. 
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Figure 41. An elevator doublet is used to excite the longitudinal dynamics of Frog. 
Test flight data is compared with simulation results to assess validity of the model. 


3. Autopilot Model Identification 

The inner-loop autopilot is a “black box” containing both sensors and con- 
troller logic. The intent was to model the unit as closely as possible without disas- 
sembling it. The published function of the autopilot is to control vertical speed and 
turn rate of the aircraft. This naturally decouples into an identification problem for 
the lateral channel, and an identification problem for the longitudinal channel. 

The lateral channel employs a rate gyro to track turn rate commands via 
feedback to the aileron. The general structure is shown in Figure 43. In order to 


identify the dynamics of the block labeled 7; in Figure 43, the unit was rotated at 


I 





Table V. Comparison of selected lateral derivatives. 


differing yaw rates. This provided a variable input signal to the feedback path with a 
frequency content covering 0 to 20 radians per second. The commanded yaw rate, r., 
was held constant. The feedback loop was broken at the summing junction of r. and 
rsp, and the feedback signal was captured. Parameter identification algorithms were 
used to determine that the feedback loop could be approximated by the following 


transfer function: 


—~.ls+ 1] 
7] = —_—_. [V.4 
s+] ( ) 


With the autopilot on, a step command in turn rate was transmitted to the 
vehicle in flight. Vehicle turn rate, as measured by the IMU, was recorded and used 
to estimate the value of Ay; at 0.25. Test data from that flight is shown in Figure 44 
and compared with simulation results using the autopilot model. 

The longitudinal channel of the autopilot senses the rate of change of static 
pressure in order to control the vehicle’s vertical velocity via feedback to the elevator. 
Flight tests data capturing the response of the vehicle to a step input in climb rate 
command was used for model identification. The identified model of the longitudinal 
channel of the autopilot is shown in Figure 45. 

With accurate models of Frog and of the onboard autopilot identified, it was 
a simple matter to determine the command bandwidths of the inner-loop controller 
(see Figure 46). For this project, this would be the limiting factor in achievable 
performance of the guidance algorithm. Hardware-in-the-loop testing of the actua- 


tors found their bandwidth to be an order of magnitude higher than the inner-loop 
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Figure 42. An aileron doublet is used to excite the lateral dynamics of Frog. Test 
flight data is compared with simulation results to asses validity of the model. 


command bandwidths. A natural next step would be to remove the autopilot after 
sufficient time and experience have reduced the risk exposure to acceptable levels, 


and exploit the extra actuator bandwidth available. 


4. Design of the Controller 


The design requirements to be met were as follows: 


1. Tracking Requirements 


e The controller should achieve perfect steady state tracking of trajectories 
in € in the presence of a constant disturbance. For a definition of the set 
€, see Chapter II. 
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Figure 43. Block diagram of the lateral channel of the inner-loop autopilot. 


2. Bandwidth Requirements 


e Command bandwidths along the lateral and longitudinal channel should be 
maximized to ensure tight tracking of the trajectory. Adequate frequency 
separation between the inner and outer loops should be assured. 


3. Closed Loop Damping and Stability Margins 


yaw rate (degrees per second} 





Figure 44. Flight test data is used to identify the dynamics of the autopilot. A step 
signal is sent to the lateral channel of the autopilot. Measured yaw rate is compared 
to simulation results. 
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Figure 45. Block diagram of the longitudinal channel of the inner-loop autopilot. 


e The dominant closed-loop eigenvalues should have a damping ration of at 
least 0.5. Phase and gain margins should be no less than 45 degrees and 6 
dB respectively in each control loop. 


4. Implementation Requirements 


e Control of Frog from the console, open-loop, must be possible in order to 
maneuver the vehicle to a suitable location overhead the field based on 
visual cues and displayed navigational data. 





Frequency [rad/sec) 


Figure 46. Bode plot of the yaw rate command to yaw rate and climb rate command 
to climb rate for Frog with the autopilot on. 
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e Switching from open-loop control to autonomous flight should be accom- 
plished from the console at the discretion of the console operator. The 
switch should be bumpless with no undesirable transients. 


e The definition of the inertial trajectory passed to the controller should be 
easily defined from the console, and expressed in terms of the parameters 
introduced in Chapter III, namely, helix angle (7.), radius (b.), and turn 


rate (tc). 


a. Linear Controller Synthesis 

The synthesis of the controller was an application of the theory outlined 
in Chapter III where the trimming trajectories were parameterized by the vector 
= lv. we ve] € R°. This was useful for the derivation of the generalized error 
dynamics. For the application at hand, however, it was convenient to remove the 
explicit dependence on v,. Since along Po € € 
DACOS(a ce 

be | 


: aT 
there is no problem in using 7, = |b. We 76| € R° to parameterize Po € E vice ne. 


— 


This was done strictly for convenience, since during the flight test, the throttle was 
to be left at a cruise power setting, and the airspeed was not explicitly controlled. 


Flight tests tracked a single trajectory in € identified as 
a 


where the units are feet, degrees and seconds, respectively. Based on (III.6), (III.7) 


and the model identification results in section 2 and 3, the nonlinear plant, 


dV = Fy(V,9,A) + Iv(V,Q)H(V,2,U), 


dt 
£Q =Fo(V,9,A) + Zo(V,Q)H(V, 2, V), 
g =} # ve ) + Lo(V, Q)H( ) (1V.6) 
ie 
ae ©), 


was available for controller synthesis. Given (IV.5) and (IV.6), the trimming values, 


Vo, Ne, and Ug were solved as the solution to 


Fy (Ve,c, Ac) + Tv(Ve, Xe )H(Ve, Qe, U.) = 9, 


ay 


Fa(Ve,Qe,Ac) + Ia(Ve,Qe JH Ve, Ve, U.) = 0, 


Q71N¢ — Ac = 0, 
bee 
ie =0 
cos(Yc) 
bepe 
- —] cos(+c) 
Co >= 0, 
B.— sin V7 


t 


ye — [0 1 0) arg(g RER) = 0, (IV.7) 





where the arg function extracts the angles X from the rotation matrix R(X): X = 


arg(R(X)). Numerically solving (IV.7), the trim values were found to be 


a 
| 


u 
| 87.7 2.037 2.43 | 


ae 
Oe | 0.00176 0.029 0.087 | 


> 
Qs 
| 


R 
| 13.41 0.023 0 , 


where, as before, standard units are feet, degrees and seconds. 
Based on lab and flight testing of the onboard sensors, the decision was 


made to design an output feedback controller of the form (see Figure IV.8): 


6H = lit 021 (oy macnn 
Oe Aan bXa + B16 kee Geos 
c ul # 1 19Acq 204-2 + Do3 (IV.8) 
£5Xc2 = 6B, 


6U = Cy6Xqy + D2b6X2 + Degd EL, 


where 6U = Ire hel. 

Classical control techniques were employed to design the controller, C;, 
shown in Figure 47. Performance objectives included setting the loop gain crossover 
frequencies at 0.5 radians per second along each channel based on the inner-loop 
command bandwidths (see Figure 46). Root-locus analysis along the lateral channel 
provides the most intuitive view of the control problem (see Figure 48). The open- 


loop plant has unstable zeros along the lateral channel due to the dihedral effect of 
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Figure 47. Linear Controller Design 


the high wing design. These unstable zeros will naturally attract the free integrators 


in the synthesis model if left unattended. 


imaginary 





Figure 48. Root-locus for the lateral channel. 


The straightforward solution was to interlace the pole-zero structure 


of the synthesis model with a stable pair of complex zeros sufficiently close to the 
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free integrators. The robustness limit was set by the dutch roll poles which were 
forced to eventually migrate to the non-minimum phase zeros. Modification of the 
dutch roll behavior would require a different approach than that chosen with the 
inclusion of the inner-loop autopilot. Replacement of the inner-loop autopilot is a 
natural extension for future flight test work. Design concerns and solutions along 
the longitudinal channel were qualitatively similar. The resulting design was a sixth 
order controller. 

Nyquist diagrams of the loop transfer function for the lateral and lon- 
gitudinal channels are shown in Figure 49 and 50. There it can be seen that phase 


and gain margin design requirements along each channel have been met. 


Nyquel Dig arn 


imaginary 





Figure 49. Nyquist diagram of the loop transfer function for the lateral channel. The 
phase margin is 112 degrees and the gain margin is 6 dB. 


5. Nonlinear Implementation 


The controller designed in the preceding section was implemented on the non- 


linear plant using the implementation proposed in Chapter III and is given next (see 
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Real 


Figure 50. Nyquist diagram of the loop transfer function for the longitudinal channel. 
The phase margin is 145 degrees and the gain margin is infinite. 


Figure 51): 
[Oy z]" = 7R(P — Po(so)), 


Eo= Ae 2 2a, 


C = Xe = AaXxe =e Bok a Bat, (IV.9) 
2X2 = CyXa + Dek + Dz E, 
U = X25: 


[be Ye We] 











Figure 51. Nonlinear Implementation of the Controller 


eh 


The following algorithm was implemented in the block labeled, Cale Mg, in 
Figure 51 in order to compute Mg. Assume the position of the vehicle is |x), y; 208 
and let the center of the trajectory be specified as [z., yc z-]’. Then, the projection 


of the position of the vehicle onto the trajectory, as defined by 77,, is 


a 
P.( so) — Ee Yo zo | ; 


where 
fy oe cos( 20%) rar. 
. , 89 COS(Yc) 
yo = db, sin(—>—~) + Yo, 
2 tebe So tan(y.) + Ze. (IV.10) 


When 7; is zero, the parameter, so, can be computed by considering the projection 


of P onto a co-level circle with radius, b., and center [z,, y, z]7. Then, 


oo (IV.11) 


Lj —L_ 


The position error resolved in the Frenet frame attached at the point Pe(so) 


is 
Tee —cos(7c) sin(#)  cos(7-) cos(#) —sin(7e) 1) — Xo 
ee = — cos( #2) — sin(;*) 0 yi — Zo | > 
Zerr sin(¥.) sin(#?) —sin(yc)cos(#)  cos(yc) Zz) — 2 
oe (1V ie 


If the helix angle is not zero, then a simple solution is not available. The point 
can be found, however, by noting that Mg must be of the form [0 ye,, Zerr|’ Then 


using (IV.12), 


— cos(Yc) sin( 5) ~ 9) + cos(7e) cos(=*)(yi — yo) — sin(ye)(21- 20) = 0, 


b, 
(IV.13) 


must hold. Using (IV.10), (1V.13) is rewritten as 


— cos(Ye) sin(*)(x) — b cos( 2oSLe)) + 2.)+ 
cos(-) cos( =2)(y1 — be sin( SPSS) 9.) 
sin(Ye)(z1 — tebeSo tania) 12.) 0, (IV.14) 


which can be solved for so. 


D. CONTROLLER IMPLEMENTATION FOR FLIGHT 
TEST 


During a flight test it 1s critical that the transition from manual to autonomous 
flight occur without any undesirable transients. The method chosen to achieve this 
was to initiate autonomous flight with the vehicle on the trajectory, and with the 
controller states recruited to their nominal values. An elegant approach toward partial 
completion of this task was to restructure the controller using the D-Implementation. 
One important convenience of this structure is the straightforward manner in which 
the controller states can be recruited at the initiation of autonomous flight based on 
the trajectory definition parameter 7,. Let t., be the time that the console operator 


turns the controller on. Then, setting 


Fo ton) = ie 
Web, tan( yc) 
0 
0 
aXe ilaee — ) (IV.15) 
0 
0 | 


perfectly recruits the controller states at the initiation of autonomous flight. The 
structure in (I[V.9) is also convenient for implementing hard limits on the feedback 


signal. It was experimentally determined that the command signals to the autopilot 
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should be kept within the following ranges 


—15 degrees per second <b, < +15 degrees per second, 


—2000 feet per minute < h < +2000 feet per minute. (IV.16) 


Since the feedback signals are the outputs of integrators in (IV.9), the states X.2 were 
easily hard limited to the values in (1V.16), with anti-windup implemented on the 
associated integrators. 

To affect a smooth transition, all that remained was to initiate autonomous 
flight when the vehicle was on the reference trajectory and aligned with it. Since there 
were no mission requirements regarding the center of the trajectory, the coordinates 
were computed to place the vehicle on the trajectory at initiation of autonomous 
flight. Furthermore, the orientation of the Frenet frame was matched to the vehicle’s 
heading. Let won be the vehicle heading at ¢,,. Let the trajectory parameter be 
specified by (IV.5). Then, setting 


cos(~) cos(7) 


Le — Lon a 9 
De 
sin(¢;) cos(Y¢ ) 
Yo = Yon — a 
Zc = Zon 


defines a trajectory where Mg at ton, as defined by (IV.12), is zero, and A3 = Ag,, 


where the subscript indicates the third element of the vector. 


1. Additional Implementation Issues 

Prior to autonomous flight of Frog, it was decided to test the RFTPS by flying 
Frog remotely from the workstation console. The console operator was provided with 
the appropriate displays showing Frog’s position, heading, and velocity, and used that 
information to command the vehicle’s turn rate and climb/descent rate in order to 
keep Frog in the local operating area. 

The vehicle’s position, as reported by the onboard GPS, is in units of latitude, 


longitude and height above mean sea level. In order to provide a more intuitive 
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navigational reference frame, the location of the workstation at the field where the 
flight tests were conducted was surveyed. This postion was used to define the origin 
of a local tangent plane. Let [Ai, \2, h]’ be the vehicle’s position, as reported by the 


GPS receiver, where 


A, = degrees of latitude, 
Ao = degrees of longitude, 
h = height GPS in meters. (1V.17) 


The variable h is a height referenced to the surface of an ellipsoid approxi- 
mating the shape of the earth. Two important parameters describing this reference 
ellipsoid, termed WGS84, are its eccentricity factor (€) and the length of its semi- 
major axis (a). The distance from the origin of the WG'S84 ellipsoid to a point on 


its surface where the local latitude is A; 1s given by 


ee (IV.18) 
1 — €?sin*(A;) 

Consider a Cartesian coordinate system with its origin at the center of the 
WGS84 ellipsoid, oriented such that its z-axis is aligned due north, its x-axis intersects 
the prime meridian, and its y-axis completes the right hand rule. This reference 
coordinate system is termed Earth Centered Earth Fixed (ECEF). Then, given the 
position reported by GPS (IV.17), the position expressed in the ECEF reference frame 


1S 


& 
| 


(N + h) cos(A2) cos(,;). 
(N + h)cos(X2) sin();). 


psa 
I 


z = (N(1—«*)+h)sin().). 


Let [\y, A2, 2]? be the surveyed position of the workstation, and let [zo, yo, zo]? be 


the same location expressed in ECEF coordinates, and suppose this position is used 


3 


to define the origin of a local tangent plane reference frame. 5ign convention and 
orientation of the local tangent plane reference frame are defined as positive x-axis 
values extending due north, positive y-axis values extending due east, and positive 
z-axis values extending straight down. Furthermore, assume [2 frog Yfrog Zfrog|. iS the 


reported position of the vehicle expressed in the ECEF coordinate system. Then, 


LI —sin(A,) cos(A;,) —sin(A2,)sin(Ai1,) cos(Aa, ) (2 frog — Com 

YI = — sin(A,) cos(,, ) 0 (frog — Yo) 

ZI — cos(Xz, ) cos(A},) + —cos(A2,)sin(Ai,) —sin(A2,) (Zfrog — Zon 
(IV.19) 


is the position of the vehicle in the local tangent plane reference frame. This position 
was displayed on the console, and later, the trajectories tracked by Frog were defined 
in this reference frame. 

The final issue that needed to be resolved was the accurate and timely calibra- 
tion of the command signals from the console to the vehicle. The autopilot onboard 
the vehicle expected to see commands in terms of PWM signals. The pulse width 
modulated signals were generated by first converting the command signals to analog 
voltages via the Digital-to-Analog converter. Then the analog voltages were directed 
to a Futaba transmitter, specially modified to respond to analog voltages instead of 
the movement of the exterior joysticks. Next, a receive module onboard Frog con- 
verted the transmitted signals to PWM format, where they were sent to the autopilot. 
A functional block diagram of the architecture used to calibrate the uplink is shown 
in Figure 52. 

The key to calibrating the uplink was the determination of the functions Fy), 
(Block 1) and Gj), (Block 2) in Figure 52. This was done in two steps. The first 
step involved sending constant commands to the DAC from the user interface (Block 
10, 11 and 12). This resulted in flight patterns consisting of numerous constant 
rate turns to the left and to the right, and to constant rate climbs and descents. 


A receive module, which was a duplicate of the one onboard Frog, converted the 
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Figure 52. Calibrating the Uplink 


transmitted signal to PWM format. This was then converted and stored in digital 
format by the RFTPS via the pulse width modulated signal converter I/O module 
(Block 3, 4 and 5). Concurrently, the onboard navigation sensors transmitted yaw 
rate and climb rate data, among other things, to the RFTPS (Block 6, 7, 8 and 9). 
Post processing of the data revealed linear relationships between the PWM signals 
sent to the autopilot and the respective climb/descent or turn rate of the vehicle 
(see Figure 53). The relationships between the command signals in terms of their 


PWM values ( and h 


and the steady-state performance of Frog in flight are 


V Coat Cpwm ) 


expressed as 


r = 0.0885r...,,, — 151.476, (IV.20) 
=: (a ca 

h = 12.3929h.,., — 18193, (IV.21) 
Saar (eee), (1V.22) 


127 


where r is the vehicle’s yaw rate in degrees per second, and h is the vehicle’s climb /descent 


rate in feet per minute. 


Yaw Rate 
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Figure 53. Top graph: measured vehicle yaw rate (degrees per second) versus PWM 
signal (4 seconds) into the lateral channel of the inner-loop autopilot. Bottom graph: 
measured vehicle climb rate (feet per minute) versus PWM signal (seconds) into the 
longitudinal channel of the inner-loop autopilot. 


The second step in calibrating the uplink involved the modified Futaba trans- 
mitter (Block 2). The known voltages applied to the transmitter were recorded, 
while the duplicate receive module was used to decode the transmitted signal, and 
the RFTPS was used to record the corresponding PWM value (Block 3, 4 and 5). 


The relationships were linear and expressed as 


Poswm = 686.7 rot, — 234.1, (IV.23) 
= Grama). 
Resum = 800Ax0% — 470, (IV.24) 
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=: ae Cc hh 


where ee is voltage sent to the longitudinal channel of the modified Futaba trans- 
mitter, and ryoi 1s the voltage sent to tlie lateral channel. 
That completed the calibration of the uplink. Given a commanded climb/descent 


rate (h.) in feet per minute, and commanded yaw rate (r.) in degrees per second at 


the console, the voltages sent to the modified transmitter were calculated as 


_ bie 
en Pa Gh Lee ie (IV .25) 
Toolt = a ie To. (IV .26) 


E. RESULTS AND ANALYSIS 


This section summarizes the results of two flight tests of the guidance and 
control algorithm developed in the preceding section. While the trajectory definition 
parameter was identical for both tests, the test flights took place on different days 
with very different environmental conditions. Furthermore, on the first test, only the 
lateral channel of the guidance algorithm was active. Longitudinal control was done 
open-loop by the console operator. Both flights took place at approximately 500 feet 
of altitude above ground level. Wind measurements were made at ground level on 
the runway, and later incorporated into the simulation. The pilot maintained control 
of the throttle throughout the tests. During autonomous flight, the pilot left the 
throttle at the trim setting for the trajectory defined. The rudder was not used. 

At the time of the tests, the RF TPS had been in operation for approximately 
three months. The two flight tests presented were not the first time the RFTPS 
was used to control Frog. Some of the model identification process, and most of 
the calibration process, used the RFTPS to control Frog. This was primarily due 
to the fact that the pilot used a conventional stick configuration to control Frog. 
Longitudinal signals corresponded to fore and aft motion, and lateral commands 
corresponding to left and right motion. Movement along one axis invariably corrupted 


the signal along the other axis. In contrast, control signals generated by the RF TPS 
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were very precise and perfectly decoupled. Also, tests had been successfully completed 
whereby a portable computer was used to send commands to the RF TPS and control 
Frog. The flight tests presented do represent the first time that a guidance and control 


algorithm was used to control Frog fully autonomously. 


1. First Flight Test 

The first test of the guidance algorithm occurred on June 20, 1997. The 
wind was measured at 10 mph out of the northwest with gusts to 15 mph. The 
pilot launched Frog and brought it up to a safe altitude. Control was passed to the 
console operator who turned the controller on. The trajectory definition 1s repeated 


for convenience, 


eee ence 
Pe 
6. = 1 1416: feeg 


) degrees per second, 


Autonomous control of the longitudinal channel was not implemented. The GPS 
operated in the differential position mode throughout the test. Prior testing of this 
GPS revealed that position error in the horizontal plane, using GPS in differential 
mode, had a standard deviation of 10 feet. The positions shown are those reported 
by GPS. The true position is unknown but thought to be within about 10 feet of that 
shown. Later, the process was repeated in simulation using the measured value of 
the wind. The trajectory flown by the vehicle is shown projected onto the horizontal 
plane in Figure 54. Additionally, the reference trajectory and the trajectory obtained 
in simulation is shown. 

The controller activity in the lateral channel is shown alongside the lateral 
error signal in Figure 55. According to the orientation of the Frenet frame, the error 
signal went negative when the vehicle was inside of the reference trajectory. The 
controller activity has noticeable high frequency content caused by the interaction 


of the one second update rate of the GPS and the slow complex zeros placed in 
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Figure 54. The first (near) autonomous flight of Frog. The commanded trajectory is 
shown as a solid line while the path flown is shown as a broken line. Also shown are the 
results from nonlinear simulation. The wind was 18 feet per second out of the west, 
northwest. The graph is oriented with the y-axis pointing due north. Approximately 
two and one half minutes of flight are shown. 


the controller. The nominal commanded yaw rate should be 5 degrees per second, 
but instead has a mean value of 8 degrees per second. This is indicative of the 
integral action compensating for modeling and calibration errors. The ground speed 
(not shown) oscillated by twice the magnitude of the measured wind. The vehicle’s 
indicated airspeed was nearly constant. Therefore, the wind appeared as a sinusoidal 
disturbance at a frequency of 0.796 radians per second. As the loop gain cross over 
frequency was 0.5 radians per second along the lateral channel, little gain was available 
to suppress this disturbance. 

Similar qualitative and quantitative performance was observed in simulation. 
The slight difference in orientation of the two trajectories is most likely due to mis- 
alignment of the modeled and actual wind. ‘Table VI quantitatively compares flight 
test data to linear and nonlinear simulation results. The close agreement in average 


(RMS) tracking error indicates that the design has good robustness properties. The 
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Figure 55. ‘Top graph: Commanded vehicle yaw rate in degrees per second from the 
controller along the trajectory, flown on June 20, 1997. Bottom graph: Lateral error 
in feet, as computed by the guidance algorithm, along the same trajectory. 


relatively large tracking errors are indicative of the low loop gain available for dis- 
turbance rejection of wind effects. In retrospect, performance requirements should 
be formulated to ensure adequate gain at these frequencies. The frequency of this 


disturbance, wg, is well known at the design stage. It can be defined using 7, as 


360 
w). COS(Y¢ ) 


Wd = 


(IV .28 


With estimates of the maximum magnitude of the wind, and requirements on the 
maximum permissible tracking errors allowed, a performance requirement could nat- 


urally be formulated as an H,, constraint. 


132 


Fs Linear Simulation | Nonlinear Simulation | Flight Test 


Lateral Mean Error -21.92 feet -22.1 feet -18.35 feet 
Lateral RMS Error 111.27 feet 111.49 feet I 24 teat 
Lateral Peak Error 140.57 feet 140.01 feet 219.54 feet 


Table VI. Error Comparison for Flight of June 20, 1997: Flight Test & Simulation 









2. Second Flight Test 

On July 23, 1997, a second flight test was conducted. The second flight test 
provided the opportunity to test the guidance algorithm in near ideal weather condi- 
tions. The wind was | to 2 miles per hour out of the west, and the air was absent of 
thermal activity and disturbances. This time, the longitudinal channel was operative. 

Frog was brought up to altitude by the pilot and control handed over to the 
console operator. The console operator switched to autonomous flight with the tra- 
jectory parameters defined as on June 20, 1997. Once again, post test flight data is 
compared with simulation results incorporating wind levels measured. The results 


are summarized in Table VII, followed by data from the flight. 


[| finear Simulation | Nonlinear Simulation | Flight Test_ 










Longitudinal Peak Error .089 feet 0584 feet ie ie etces 


Table VII. Error Comparison for Flight of July 23, 1997: Flight Test & Simulation 





The resulting path in space flown on July 23 is shown in Figure 57. It is 
compared to the reference trajectory. The projection of the trajectory tracked onto 


the horizontal plane is shown in Figure 56. It, also, is compared to results from 
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simulation. The history of the controller activity is shown for the lateral channel in 


Figure 58, and for the longitudinal channel in Figure 59. 





Figure 56. The second autonomous flight of Frog. Weather conditions were near 
ideal as Frog tracked the 3-D trajectory shown by the dashed line. The results from 
simulation are shown as a dash-dot line. 





Figure 57. The second autonomous flight of Frog. Weather conditions were near ideal 
as Frog tracked the 3-D trajectory shown. 
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Figure 58. Top graph: Commanded vehicle yaw rate, in degrees per second from the 
controller, along the trajectory flown on July 23, 1997. Bottom graph: Lateral error 
in feet, as computed by the guidance algorithm, along the same trajectory. 


With little wind, the test data characterizes the performance of the guidance 
and control algorithm in the presence of modeling errors, transport lag, hardware non- 
linearities, and sensor noise. The simulation results were a close match to flight test 
data along the lateral channel. The average tracking error just over 21 feet. It can be 
seen in Figure 58 that the integral control is holding about one half degree per second 
commanded yaw rate to compensate for constant disturbances. In flight, tracking 
errors along the longitudinal channel were even less than lateral tracking errors. The 
average command signal along the longitudinal channel was —5 feet per second. This 
probably indicates that the power setting was too high. Subsequent testing will 
incorporate airspeed and throttle control which will address this issue. In simulation, 
however, longitudinal tracking errors were extremely small. The discrepancy is most 
likely the result of unmodeled vertical disturbances in the airmass due to light thermal 


activity. 
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zettor (feel) and commanded cimb rate {Ips) 





Figure 59. The output of the longitudinal channel of the controller is shown along 
with the error signal on which it is acting. Note that the positive z-axis points down 
and that positive error corresponds to the vehicle being below the reference trajectory. 


F. CONCLUSIONS 


In this chapter, an integrated system for the design, development and testing of 
guidance, navigation, and control algorithms for unmanned air vehicles was presented. 
Extensive use was made of commercial off-the-shelf (COTS) hardware. This kept costs 
low, and made the system easily scaleable. Sophisticated autocode tools allowed a 
two man team to write, test, and maintain thousands of lines of error free real time 
code. In a single, unified environment, the avionics system was designed, simulated, 
simulated incorporating hardware-in-the-loop (HITL), tested in flight, and used to 
control an aircraft in flight. 

As a proof-of-concept demonstration, the RFTPS was used to take the in- 
tegrated guidance and control concept presented in Chapter III and evaluate it in 
the environment for which it was proposed to be used. The project began with an 
unmanned air vehicle with unknown flight characteristics. In addition, a “black box” 
autopilot with unknown internal dynamics was placed onboard the vehicle. Through 


the use of the RF TPS, a high fidelity simulation of the vehicle and autopilot was built 
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and used to synthesize the applicable guidance and control laws. Extensive testing 
in simulation followed, and appropriate user interfaces were developed. Then, the 
RFTPS was used to control the vehicle in flight using the guidance and control laws 
developed, collect the data during the flight test, post process the data, and evaluate 
the performance of the algorithms. 

The success of the project demonstrated both the utility of the integrated guid- 
ance and control algorithm as well as the capabilities of the RFTPS. The integrated 
guidance and control algorithm was shown to work well in a real world application. 
Performance in flight was very close to performance in simulation, which speaks well 
of the robustness properties of the controller. Changing environmental conditions 
highlighted some disturbance rejection issues that should be addressed in the future. 
The RFTPS was shown to be powerful, portable, rugged, effective in the field and in 


the lab, and safe and reliable at controlling an aircraft. 
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APPENDIX A. THE TAIL SIZING DESIGN 
TOOL 


1. INTRODUCTION 
This appendix details the operation of the Taz Sizing Design Tool. The Tail 


Sizing Design Tool runs from the MATLAB command line, and requires the LMI and 
CONTROL Toolboxes. The purpose of the Tail Sizing Design Tool is to allow the 
user to map and view the Tail Sizing Design Space for a given aircraft definition, or 
compare the Tail Sizing Design Space of two different aircraft definitions. There are 
three versions of the Tail Sizing Design Tool. The Tail Sizing Design Tool Version 
A maps the Tail Sizing Design Space for the high angle-of-attack excursion search- 
ing over static, state feedback, controllers. Version B solves the same problem but 
searches over dynamic, output feedback, controllers. Version C maps the Tail Sizing 
Design Space for the wind-shear penetration constraint searching over static state 
feedback controllers. Version A and B are well developed and allow the user to in- 
clude aeroelastic effects, the addition of a canard, and to view the Tail Sizing Design 
Space in a three dimensional plot. Version C is less well developed and provides for 
only rigid body, tail only, aircraft dynamics, and two dimensional views. 

The general structure of all three versions is similar. The Tail Sizing Design 
Tool is comprised of three directories which must reside on the same level within 
the user’s file system. The first directory is termed GUI, which is an acaronym for 
graphical user interface. As the name suggests, it contains files which allow the 
user to specify an aircraft definition, generate a database for an aircraft definition or 
view a Tail Sizing Design Space for an aircraft definition by pushing the appropriate 
“buttons” on the display generated. The second directory is termed MODEL, and 
it contains the files necessary to model the dynamics of an aircraft. This directory 
contains the user supplied aerodynamic stability and control derivatives. The third 


directory is termed DATA, and it contains data files generated by running the Tail 
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Sizing Design Tool. 
Assume for the moment that the proper stability and control derivatives have 


been placed appropriately in the MODEL directory. The process is as follows. 


1. From the GUI directory, run s2ze_it. 


i) 


. Fill in the appropriate values for the range of tail volumes and actuator rates 
of interest. 


3. Specify the aircraft definition parameters such as the use of a canard or inclu- 
sion of a flexible mode. 


4. For output feedback controllers, specify the output sensors to be used from 
the available choices. 


. Name the database file and start the generation of the data. This executes 
the program getdata. 


ot 


6. Run the program view_it to querry a database previously generated. 


2. GENERATING A DATABASE 


First, we will consider the two programs szze_zt and getdata. The interface 
generated by size_it (Version B) is shown in Figure 60. By selecting the push button, 
"Output Page,” another GUI is created (Figure 61). This provides the user the option 
of specifying which sensors to use for feedback. When the user is satisfied with the 
definition of the aircraft, he begins execution by selecting the pushbutton Run Itl. 
This executes the program Gretdata, shown next for Version B. 
hhihhhhhhhhhhhhhhhhhhhhhhhhhhhh 
if Getdata.m i 
WNhhhhhhhhhhhhhhhhhhhhhihhhhhh 


4 Initialize global variables and move to the model directory. 
cgrate=[]; data=[]; 

oie ener 

cd model 


Use the graphics handles to transfer the user defined aircraft 
Adefintions to variables in the workspace. Begin with the frequency 
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Figure 60. Interface called by sizeit.m 


4% of the flexible mode. 
freqi = str2num(get (freq_ed,’string’)); 
save flexdata freql 
4This is the location of the flexible mode sensors. 
sens_loc = str2num(get(sens_ed,’string’)); 
4Some locations (nodes/anti-nodes) make the problem 
Aunobservable. Check for this and alert the user if it is a problem. 
[phi , phidot ]=modeshpe(sens_loc4250/100) ; 
if (phidot < te-5), 
disp(’Flexible pitch rate is nearly unobservable. ’) 
disp(’Change the sensor location ... stopping’) 
elf 
/, Message 
frame_title = uicontrol(gcf,... 
estyle- text’ ,... 
Bpositaone 120.270 200 30],.... 
*string’,’There was a problem,’,... 
mecweproundcolor’,[11 1],... 
wpackeroundcollor’ , [0 0 0] ); 
4 Message 
frame_title = uicontrol(gcf,... 
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Figure 61. Selection of output sensors 


SUV Les = Ge xan gener 
’position’ , [120 240 200 30], 
*string’,’see the command window.’,... 
’foregroundcolor’,[1 1 1],... 
’backgroundcolor’,[0 0 0]); 

clear 

break; 
end; 
Name the file to store the data. 
opfname = get(save_ed,’string’) ; 
AInteger 1 if model is aeroelastic. 
flexflag = get(ck_flex,’value’); 
AInteger 1 if the model has a canard. 
canardflag = get (ck_canard, ’value’) ; 
Alf a canard is present, specify its canard volume. 
if (canardflag==1), 
split = str2num(get (canard_ed,’string’)); 
else, split = 0; end; 
“This is the range and increments of tail volume to check. 
vhlo = str2num(get(minvol_ed,’string’ )) ; 
vhhi = str2num(get(maxvol_ed,’string’)) ; 
vhinc = str2num(get (incvol_ed,’string’)); 
“This is the required closed loop damping. 
damping = str2num(get (damp_ed, ’string’)); 
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theta = 2%acos(damping); % pole placement region definition 
“This is the range and increment of actuator rates. 
ratelo = str2num(get(minrate_ed,’string’)) ; 

ratehi = str2num(get (maxrate_ed,’string’)); 

rateinc = str2num(get(incrate_ed,’string’)); 

%AThis specifies which sensors to use for output feedback. 
4, An integer 1 indicates use of the sensor. Choices are 
74 flight path angle, airspeed angle of attack, local 

4, pitch rate, local pitch angle. The default is full 

7, information output feedback. 

gamaflag = get(ck_gama, ’value’) ; 

asflag = get(ck_as,’value’); 

aoaflag = get(ck_aoa,’value’); 

lprflag = get(ck_lpr, ’value’) ; 

lpaflag = get(ck_lpa,’value’) ; 

fsflag = get(ck_fs,’value’) ; 


4, The actuator dynamics are specified in the function 
7% ACTAUTOR.M. Canard and tail actuators do not have 

4, to be the same although dummy dynamics must be 

4 specified for a canard even if not used. 
[Aa,Ba,Ca,Da] = actuator; 

act=length(Aa) ; 

4, For the first order actuator modelused here, 

/4, specify the time constant 

a_bw = 20; 


4 The angle of attack excursion is equal to acos(gust/vt). 
gust = 66; fps 


% The equilibrium point is specified by the flight 
YZ, path angle and true airspeed of the HSCT. 
gamma = 0; vt = 422; 4% deg fps 


4 Specify tolerance limit for the cg binary search. 
tol = 0.01; 


4 MATLAB uses the matrix variable ’’region’’ to define 
4 the closed loop poleplacement region. Form it 

4 based on user inputs. 

region=zeros (4,8) ; 

region(1,1)=.02+11; 
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region(2,2)=-40+11; 
region(3,3)=0+2i; 
region(1,5)=1; 
region(2,6)=-1; 
region(3,7)=cos(damping) ; 
region(4,8)=cos(damping) ; 
region(3,8)=-sin(damping) ; 
region(4,7)=sin(damping) ; 


74 Set up three loops to map the design space. 
/, First, loop through range of tail volumes. 
for vh=vhlo:vhinc:vhhi; 


7, Initialize some variables for each tail volume. 
yemax=0; ydemax=0; ycmax=0; ydcmax=0; 


/, Loop through the range of actuator rates. 
for ratelimit = ratelo:rateinc:ratehi; 


7, Reset cg upper and lower bounds for the 
4 binary search. 
xcgmax = 1.0; xcgmin = 0; 


4 Binary search on c.g. ;drop out when 
4 within tolerance. 
while (norm(xcgmax-xcgmin) > tol); 


save cgvh xcg vh split; 


4 Load appropriate data for trim; initial 
4% guesses for trim. 
load trimdata 


4 Trim at cg specified and linearize at peak 

4 of gust condition. 

if ( flexflag == 1), 

[xtrim,utrim,ytrim] = trim(’eom_b70w’ ,x0,u0,.... 
[vt;gamma],(],(1],{1;2]); 

xtrim(2) = xtrim(2) + gust; 

laf ,bf,cf,df] = linmod(’eom_b70w’ ,xtrim,utrim) ; 
bf = biC2 12 2e die — a 
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mecamee) = xtrim(2) = gust; 

else, 

rnin ucrim,ytrim) = trim(’eom_b70r’ ,xOr,u0,.... 
moeeanmal, fi), (1), 01;21); 

xtrim(2) = xtrim(2) + gust; 

femebt cr dt) = linmod(’ecom_b/0r’ ,xtrim,utrin) ; 
Dae br: 1:2); df = df(: ,1:2); 

mebim 2) = xtrim(2) - gust; 

end; 


Waiconmect the actuators in series. 
[A,B,C,D] = series(Aa,Ba,Ca,Da,af,bf,cf,df) ; 
[m,n]=size(A); 


%, Place the actuator states after the plant states. 
(n,p]=size(B) ; 

T1 = [zeros(n-2,2) eye(n-2);eye(2) zeros(2,n-2)]; 
eee 1 ,A/T1; 

B = T1%B; 

CUT i- 


©2 
ul 


4 Similarity transformation used to make local 
4 pitch rate and local body angle states vice 
4% using generalized coordinates. 
if (flexflag==1), 
[phi, phidot ]=modeshpe(sens_loc%250/100) ; 
T2 = Leye(4) zeros(4, 2+act) ; 
(0 0 0 1 phidot 0] zeros(1,act); 
[ORG eOeGaphidot] zeros (1, act); 
zeros(act,6) eye(act)]; 
else, T2=eye(n); end; 


Abar = T2A/T2; 
Bbar = T2B; 
Cbar = C/12- 


4 Form output matrix based on user specified sensors. 
4% airspeed is Cbar(1,:) 

4 flight path angle is Cbar(2,:) 

%Z angle of attack is Cbar(3,:) 

4 local pitch rate is 5th state [0000100 0] 
4 local pitch angle is 6th state [0000010 0] 
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Cuseclleer— ale 

if gamaflag == 1; 

C_ select = "(6_ select Chamero nn 
end; 

if asflag == 1; 

C_select = [Clselleet, CboanGuaa tn 
end ; 

1f aoaflag == 1; 

G_select = (Clselect. Charis on 
end; 

1f lprflag == 1; 

C_select = [C_select;[zeros(1,n-4) 10 0 0]]- 
end; 

1f lpaflag == 1; 

C_select = [C_select wizeros@in-2e0 175070] 5 
end; 

if fsflag == 1; 

C_select = eye(n); 

end; 
Cbar=C_select; 

[n,m]=size(Bbar) ; 
[p,n]=size(Cbar) ; 


4, Use the angle of attack excursion to define 
7, the initial condition. 

i0=[0 gust 0 0 zeros(1,n-4-act) zeros(1,act)]’; 
4i0=T%i0; 


4 Specify rate and ampltude constraints. The 
7, rate limit comes from the 
4 GUI. The saturation limit is set at 30 degrees. 


ucmax = 30; uemax = 30; 
udcmax = ratelimit; udemax = ratelimit 


4 Finally, form the LMIs. 
setlmis([]); 


R 
5 


imivar(1, Pm 1]) 39) Ry 
imivar(i, [in 4) oe, 
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iimvvar GZ meinen); wAeAk 7 
ead Zemin pl); wy, Bk % 
imivar (2. an mn) ); ole i, 


4 This is pole placement LMIs formed using 
4 LMI region matrix REGION. 
[rr,rc]=size(region) ; 

iene =241T, 

error(’REGION must be a Mx(2M) matrix’); 


end 
ri=0; 


4 scan each region (L,M) (diagonal blocks) 
while ri<rr, 
bs=round(imag(region(rit+1,ri+1))); 


fe bs, bs=rr-ri; end 


ereal(rerion(riti:ritbs,ri+1:ri+bs)) ; 
M=region(ri+1:ritbs,rrtrit+1:rrt+ritbs) ; 
pole=newlmi; % pole placement in the region (L,M) 
for 1=1:bs, 


nbi=2/4i-1; % row of block in target LMI 
4, off-diagonal 2x2 block 
Bor joi: 1-1; 


nbj=2/4j-1; % col of block in target LMI 
lmiterm([pole,nbi,nbj,R] ,L(i,j),1); 
lmiterm([pole,nbi,nbj,R] ,Abar,M(i,j)); 
Imaterm( [pole,nbi,nb;,R) ,M(j,1),Abar’); 
lmiterm([pole,nbi,nbj,Ck] ,Bbar,M(i,j)); 
lmiterm([pole,nbi,nbj,-Ck] ,M(j,i),Bbar’) ; 
lmiterm([pole,nbit1,nbj+1,S] ,L(i,j),1); 
lmiterm([pole,nbi+1,nbj+1,S],M(i,j) ,Abar) ; 
lmiterm([pole,nbi+1,nbj+1,S],Abar’,M(j,i)); 
lmiterm([pole,nbit1,nbj+1,Bk] ,M(i,j) ,Cbar) ; 
lmiterm([pole,nbit1,nbj+1,-Bk] ,Cbar’ ,M(j,i)); 
lmiterm([pole,nbiti1,nbj,0],L(i,j)); 
lmiterm([pole,nbi+1,nbj,Ak] ,M(i,j),1); 
lmiterm([pole,nbit+1,nbj,0],M(j,i)%Abar’); 
lmiterm([pole,nbi,nbj+1,0],L(i,j)); 
lmiterm( [pole ,nbi,nbj+1,-Ak] ,M(j,i),1); 
lmiterm([pole,nbi,nbj+1,0],M(i,j)%Abar) ; 


end 
7, diagonal 2x2 block 
lmiterm([pole,nbi,nbi,R] ,L(i,i),1); 
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lmiterm([pole,nbi,nbi,R] ,Abar,M(i,i),’s’) 
lmiterm([pole,nbi,nbi,Ck] ,Bbar,M(i,i),’s’ 
lmiterm( [pole ,nbi+1 ,nbi+1,S],L(i,i),1); 
lmiterm([pole,nbit1,nbit+1,S],M(i,i),Abar,’s’); 
lmiterm([pole,nbit+1,nbi+1,Bk],M(i,i),Cbar,’s’); 
lmiterm([pole,nbit+1,nbi,0],L(i,i)); 
lmiterm([pole,nbit+1,nbi,Ak] ,M(i,i),1); 
lmiterm([pole,nbi+1,nbi,0] ,M(i,i)/Abar’); 

end 

ri=ritbs; 

end % while 


« 

? 
ye 
3 


%, This is x_0 LMI from the invariat elipsoid. 
xO_lmi = newlmi; 

lmiterm ( [-xO_lmi 1 1 0], 1); 

lmiterm( [-xO_lmi 2 1 

lmiterm( [-x0O_lmi 2 2 

Imiterm( [—x0limi Sede oil eo 
lmiterm( [-xO_lmi 3 2 0], eye(n)); 
Imi tern G@eatex0elma SS 5) aoa 


7, This is (u_max) saturation limit elevator LMI. 
umaxe_lmi = newlml1; 

lmiterm([-umaxe_lmi 1 1 R], 1,1); 

lmiterm( [-umaxe_lmi 2 1 0],eye(n)); 
lmiterm([-umaxe_lmi 2 2 S], 1,1); 
lmiterm([-umaxe_lmi 3 1 Ck], [0 1i],1); 
lmiterm([-umaxe_lmi 3 2 0], Of*eye(n)); 

lmiterm( [-umaxe_lmi 3 3 0], uemax2); 


4, This is (u_max) saturation limit canard LMI. 
if (canardflag==1), 

umaxc_lmi = newlmi; 

lmiterm([-umaxc_lmi 1 1 R], 1,1); 

lmiterm( [-umaxc_lmi 2 1 0],eye(n)); 
lmiterm([-umaxc_lmi 2 2 S], 1,1); 
lmiterm([-umaxc_lmi 3 1 Ck], [1 0],1); 
lmiterm([-umaxc_lmi 3 2 0], O%eye(n)); 
lmiterm( [-umaxc_lmi 3 3 0], ucmax*2); 

end; 
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%, This is (ud_max) actuator rate elevator LMI. 
udmaxe_lmi = newlmi; 

Imiterm([-udmaxe_lmi i i RJ, i,i); 

lmiterm( [-udmaxe_lmi 2 1 0],eye(n)); 
Imiterm([-udmaxe_lmi 2 2 S], 1,1); 
lmiterm([-udmaxe_lmi 3 i Ck], [0 a_bw],1); 
Imiterm([-udmaxe_lmi 3 i R], [zeros(i,n-1i) -a_bw],1); 
Imiterm([-udmaxe_lmi 3 2 0], [zeros(i,n-1) -a_bw]); 
lmiterm( [-udmaxe_lmi 3 3 0], udemax~2); 


fetnis as (ud _max) actuator rate canard LMI. 

if (canardflag==1), 

udmaxc_lm1 = newlmi; 

Imiterm([-udmaxc_lmi i i RJ], 1,1); 

Imiterm( [-udmaxc_lmi 2 1 0],eye(n)); 
Imiterm([-udmaxc_lmi 2 2 S], 1,1); 
Imiterm([-udmaxc_lmi 3 i Ck], [a_bw 0],1); 
lmiterm([-udmaxc_1mi 3 1 R], [zeros(1i,n-2) -a_bw 0],1); 
Imiterm([-udmaxc_lmi 3 2 0], [zeros(i,n-2) -a_bw 0]); 
Imiterm( [-udmaxc_lmi 3 3 0], udcmax~2); 

end; 


2 
1 
1 
Z 


lmisys = getlmis; 


pecwech to see if the set of LMIs is feasible. 
4 TMIN should be negative for feasibility and 
4, XFEAS is the decision vector that validates 
7, the above LMIs. 

[tmin,xfeas] = feasp(lmisys, [0 200 ie8 0 0]); 


4 Values of TMIN less than 0 are not strictly feasible 
% but seem to be close enough to work. 
if tmin < ie-4; 


7, Recover the LMI variables. 
R = dec2mat (lmisys,xfeas,R) ; 
S = dec2mat (lmisys,xfeas,S) ; 


ak = dec2mat (lmisys,xfeas,Ak) ; 
bk = dec2mat (lmisys ,xfeas ,Bk) ; 
ck = dec2mat (1lmisys ,xfeas ,Ck) ; 
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% Construct the controller and form the close loop system. 
[u,sd,v]l=svd(eye(n)- R%4S); % factorize I-RS 
sd=diag (sqrt (1./diag(sd))); 

Ni=sd*v’; Mti=u*sd; 

Ac=Nix* (ak-S* (Abar) R-bkCbarR-S*Bbar*ck) *Mti; 
Bc=Nix (bk) ; 

Cc=(ck) *Mti; 

P = ltisys(Abar, Bbar,Cbar,zeros(p,m)) ; 

cont = ltisys(Ac,Bc,Cc,zeros(m,p)); 

clsys = slft(P,cont,m,p); 

[acl,bel,c¢cel,dcll] = itiss(elsysye 


4 Simulate the closed loop system with the specified 

% initial condition and record the maximum actuator 

7, amplitudes and rates. 

[ye,x,t] = initial (acl, [Bbar(: ,2) ;zeros(@i. eee 
[zeros(1 ,n) [0 1) *Cel- 0 (0 3zer6s (ns Dee 

yemax = max(abs(ye)) ; 

Lyc,x,t] = initial (acl, [Bbar(:.2) > zerostn. lee 
[zeros(?,n) [17 0]*Cel{0 0: Zeros Ge) pr 

ycmax = max(abs(yc)); 

lyde,x,t] = initial (acl; Bbar(- 42) .zeroc i eee 
[[zeros(1,n-1) -a_bw] [0 a_bw]*Cc] ,0,[10;zeros(n,1)]); 
ydemax = max(abs(yde)) ; 

[ydc,x,t] = initial (acl, [Bbar( 42) zeros. 1) 
[[zeros(1,n-2) -a_bw 0] [a_bw 0]*Cc] ,0,[10;zeros(n,1)]); 
ydcmax = max(abs(ydc)); 

eigvalue = eig(acl) 

data = [data; xcg ratelimit vh split ycmax yemax ydcmax.... 
ydemax tmin utrim’ xtrim’ eigvalue’]; 


4 Move cg limit aft. 
xcgmin = xcg; 


4 Otherwise, move cg limit forward, 
else, 

xcgmax = xcg; 

end; ASET FEASIBLE, if 


% Chech to see if aft cg limit is 


4 within tolerance 
if ( abs(xcgmax-xcgmin) < tol | xcg > .95 | xcg < .05 ) 
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break; 
end; Y%WITHIN TOLERANCE, if 
end; %ZBINARY CG SEARCH, while 


7, This is the aft most cg for the rate, tail 

7% volume combination. 

if (ydcmax ~= 0 | ydemax ~= 0), 

cgrate=[cgrate; vh xcg ydcmax ycmax ydemax yemax eigvalue’]; 
end; 

end; ARATE SWEEP, for 

end; %VH SWEEP, for 


/, Change to the DATA directory and store the 
7, result with the name provided from the GUI. 
4, Return the user to the GUI directory. 
ea 
cd data 
eval([’save ’,opfname,’.mat cgrate data ’]) 
eG 2.3 
cd gui 
The above code is associated with the output feedback problem. For the state 
feedback problem, the code is very similar. The essential difference is in the LMI 


formulation and reconstructionof the controller. That part of the code from Version 


A substantially different from Version B is shown next. 
7, Form the state feedback LMIs 


setlmis([]); 


~< 
i 


selmivar(t. in 1]): 
iimaiais@2 sho nil): 


= 
i" 


4 This is Lyapunov lmi with real part constraint 
lyap = newlmi; 

lmiterm([lyap 1 1 Y], Abar+tcenter*eye(n), 1, ’s’); 
lmiterm([lyap 1 1 W], Bbar, 1, ’s’); 


4 This is the pole placement LMI 


15] 


pplac = newlmi; 


lmiterm([pplac 1 1 YJ], sin(theta/2)*Abar, 1, ’s’); 
lmiterm({pplac 1 1 W], sin(theta/2)*Bbar, M, ’s’); 
lmiterm({[pplac 2 1 Y], cos(theta/2)*Abar, 1); 
lmiterm([pplac 2 1 W], cos(theta/2)*Bbar, M); 
lmiterm([pplac 2 1 -Y], -cos(theta/2),Abar’) ; 
lmiterm({[pplac 2 1 -W], -cos(theta/2)*M’ ,Bbar’) ; 
Imiterm([pplac 2 2 Y], sin(theta/2)*Abar, 1, ’s’); 
lmiterm({[pplac 2 2 W], sin(theta/2)*Bbar, M, ’s’); 


% This is the Y > O LMI 


ypos = newlmi; 
Imiterm({-ypos 1 ig ies 


7, This is x_0O LMI 

xO_lmi = newlmi; 
Imiterm({-xO_lmi 1 1 0], 1); 
Imiterm( l=xO02lmi 26160 e10)< 
Imiterm([-x0_lmi 2 2 Y],1,1); 


4, This is ve_max LMI 


uemax_]lmi = newlmi; 

lmiterm({-uemax_lmi 11 Y], 1,1); 
lmiterm({-uemax_lmi 2 1 W], C_amp_e,M); 
lmiterm([-uemax_lmi 2 2 0], (uemax2)); 


%, This is ude_max LMI 


udemax_lmi = newlmi; 


lmiterm([-udemax_lmi 1 1 YJ], eye(n),eye(n)); 
lmiterm([-udemax_lmi 2 1 W], C_rate_e*Bbar,M); 
lmiterm({-udemax_lmi 2 1 Y], C_rate_e*Abar,eye(n)); 
lmiterm({-udemax_lmi 2 2 0], udemax~2); 


if (canardflag == 1), 
4 This is uc_max LMI 


ucmax_lmi = newlmi; 
Imiterm([-ucmax_lmi 11 Y], 1,1); 


lmiterm([-ucmax_lmi 2 1 W], C_amp_c,M); 
lmiterm([-ucmax_lmi 2 2 0], (ucmax~2)); 


7, This is udc_max LMI 

udcmax_Ilmi = newlmi; 

lmiterm([-udcmax_Imi 1 1 Y], eye(n),eye(n)); 
lmiterm([-udcmax_lmi 2 1 W], C_rate_c*Bbar,M); 
Imiterm([-udcmax_lmi 2 1 Y], C_rate_c*Abar,eye(n)); 
Imiterm([-udcmax_lmi 2 2 0], udcmax~2); 


end; 


lmisys = getlmis; 


[tmin,xfeas] = feasp(lmisys,[0 200 1e9 0 0]); 
1f tmin < le-4 


4, recover LMI variables 


y = dec2mat (lmisys,xfeas,Y) ; 
w= dec2mat (lmisys,xfeas,W) ; 
Kbar = w*M/y; 


% record initial condition responses 


if (canardflag == 1), 

yc = initial (Abar+Bbar*Kbar,Bbar,C_amp_c*.... 
Kper 00 0] ,10); 

ycmax = max(abs(yc)); 

yre = initial (Abar+Bbar*Kbar,Bbar(:,1),... 
C_rate_c*[Bbar*Kbar + Abar],0,i0); 

ydcmax = max(abs(yrc)) ; 

end; 


ye = initial (Abar+Bbar*Kbar , Bbar,C_amp_e*Kbar, [0 0] ,i0); 
yemax = max(abs(ye)); 

yre = initial (Abar+Bbar*Kbar ,Bbar(:,2),C_rate_e*.... 
[Bbar*Kbar + Abar],0,i0); 

ydemax = max(abs(yre)) ; 

elgvalue = eig(Abar+Bbar*Kbar) ; 
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data = [data; xcg ratelimit vh split ycmax yemax .. 
ydcmax ydemax tmin utrim’ xtrim’ Kbar(1,:) Kbar(2,:) eigvalue’]; 


Again, for the gust penetration/state feedback problem (Version C), the code 
is very similar. Principally, the formation of the lmear model (which includes gust 
dynamics) is different. After that, the code proceeds in a similar manner to the angle 
of attack/state feedback problem (Version A). That part of the code from Version C 


substantially different from Version A is shown next. 


7%, Gust penetration model formulation, 
% including gust dynamics and 
7%, LMI formulation. 


% Load appropriate data for trim 

load trimdata 

([xtrim,utrim,ytrim] = trim(’eom_b70r’,x0r,u0,.... 
[vt;gamma],[],[1],(1;2]); 


/, Linear model for given cg/vh 


laf ,bf,cf,df] = linmod(’eom_b70r’ ,xtrim,utrim) ; 
jones one AS ae = Che 3 2) 


7, Connect the actUators in series 
[A,B,C,D] = series(Aa,Ba,Ca,Dasaf. bi ct dt)- 
c=[]; D=(]; 


% Place the actuator states after the plant states 
[n,p]=size(B); 
T = [zeros(n-1,1) eye(n-1);eye(1) zeros(1,n-1)]; 


A = T*A/T; 
B = TB; 
C=eye(n); 


D=zeros(n,p); 


7, Add gust states to end of model 


Lo4 


Ag=zeros(n+2 ,n+2) ; 

Ag(n+1,n+1)=lambdai ; 

Ag (n+2,n+2)=lambda2 ; 

fel. 5,nt1)=-A(1:5, 2); 

eagle Oo n+2)=-A(1:5, 2) ; 

Abari=[A zeros(n, 2) ;zeros(2,n+2)] + Ag; 
Bbar1=(B;0;0]; 

Cbar=eye(nt2) ; 

Dbar=zeros(n+2,p) ; 

T = eye(n+2) ; 

T(2,nt+1) = -1; 

T(2,n+2) = -1; 

Abar=T*Abar1/T: 

Bbar=T*Bbar1; 

Pespecity Static Controller Structure 

M = [eye(n) zeros(n,2); zeros(2,n+2)]; 

4 Define gust i.c. 

i0=[0 0 0 0 zeros(1,n-4-act) zeros(1,act) -gust gust]’; 
7 Rate/Amplitude selection matrices 
C_amp_e = 1; C_rate_e = [zeros(1,n-act) 1]; 
4 modify for gust states 

(erate ec — [CG rate_e 0 0]: 

7 AOA output matrix 

C_alpha = (57.3/vt)*(0 1 zeros(1,n-2) 0 0]; 
4 Specify rate and amplitude constraints 
uemax = amplimit; 


udemax = ratelimit; 
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% form the LMIs 
setlmis([]); 


fi 
W 


Imivar(i,{n 1;2 1]); 
Imivar(2,[1 n+2]); 


=~ 


This is Lyapunov 1lmi with real part constraint 
lyap = newlmi; 

lmiterm([lyap 1 1 Y], Abart+tcenter*eye(n+2), 1, ’s’); 
lmiterm([lyap 1 1 W], Bbar, M, ’s’); 

%Z This is the pole placement LMI 


pplac = newlmi; 


lmiterm([pplac 1 1 Y], sin(theta/2)*Abar, 1, ’s’); 
lmiterm([pplac 1 1 WJ], sin(theta/2)*Bbar, M, ’s’); 
lmiterm([pplac 2 1 Y], cos(theta/2)*Abar, 1); 
lmiterm([pplac 2 1 W], cos(theta/2)*Bbar, M); 
lmiterm([pplac 2 1 -Y], -cos(theta/2) ,Abar’) ; 
lmiterm([pplac 2 1 -W], -cos(theta/2)*M’,Bbar’); 
lmiterm([pplac 2 2 Y], sin(theta/2)*Abar, 1, ’s’); 
lmiterm([pplac 2 2 W], sin(theta/2)*Bbar, M, ’s’); 


Y This is the Y >40) EMI 


ypos = newlmi; 
Imiterm([-ypos 1 1 ¥],1,1); 


io Thiseis x20 oL ME 

xO_lmi = newlmi; 
lmiterm([-x0O_1lmi 1 1 0], 1); 
lmiterm([-xO_lmi 2 1 0], i0); 
lmiterm([-xO_lmi 2 2 Y],1,1); 
%4, This is ue_max LMI 
uemax_lmi = newlmi; 


Imiterm([-uemax_lmi 1 1 Y], 1,1); 
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Imiterm([-uemax_lmi 2 1 W], C_amp_e,M); 
lmiterm([-uemax_lmi 2 2 0], (uemax~2)); 


% This is alpha_max LMI 


Jamax_lmi = newlmi; 

%Zlmiterm([-amax_lmi 1 1 Y], 1,1); 
Zlmiterm([-amax_lmi 2 1 Y], C_alpha,1); 
Zlmiterm((-amax_lmi 2 2 0], (alphalimit~2)); 


7, This is ude_max LMI 


udemax_lmi = newlmi; 

lmiterm([-udemax_lmi 1 1 Y], eye(n+2) ,eye(n+2)); 
lmiterm([-udemax_lmi 2 1 W], C_rate_e*Bbar,M) ; 
lmiterm([-udemax_lmi 2 1 Y], C_rate_e*Abar,eye(n+2)); 
Imiterm([(-udemax_lmi 2 2 0], udemax~2); 


lmisys = getlmis; 

[tmin,xfeas] = feasp(lmisys,[0 200 1e9 0 0]); 
if tmin < ie-4 

“4% recover LMI variables 

y = dec2mat (lmisys,xfeas,Y) ; 

w = dec2mat (lmisys,xfeas,W) ; 

Kbar = w*M/y; 

% record initial condition responses 
ye = initial (Abar+Bbar*Kbar,Bbar,.... 
C_amp_e*Kbar ,0,10) ; 

yemax = max(abs(ye)) ; 

yre = initial (Abar+Bbar*Kbar,Bbar,.... 
C_rate_e*[Bbar*Kbar + Abar] ,0,i0); 


ydemax = max(abs(yre)) ; 


eigvalue = eig(Abar+Bbar*Kbar) ; 
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When the process has completed for the range of values specified for either Version 
A, B, or C, the user is returned to the MATLAB command prompt. The program 
results will have been stored in the DATA directory using the name supplied by the 


user. The user is free to view them or generate another data file. 


3. VIEWING A DATABASE 


After a database has been generated for a particular aircraft definiton, the 
program view_it.m is used to view the data. The program view_it.m presents the 


graphical user interface shown in Figure 62. 
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Figure 62. Interface called by viewit.m 


The interface is fairly self explanatory. The options are to view the three 
dimensional Tail Sizing Design Space, or to view a 2D slice of the Tail Sizing Design 


space. Additionally, the user can specify a second database, in order to compare two 
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different aircraft definitions on a single graph. If a two dimensional view is desired, 
then the user must specify the actuator rate limit used to intersect the lower surface 
of the Tail Sizing Design Space. When ready, the user selects the button labeled 


” Show It,” which executes the program in the GU/ directory called showdata.m. 


i hihkhbhihhhhbbhhhhhhbhheh 
7, Showdata.m i, 
oh ioAiihhittetotrhtttel htt hiehtttete 
yi=[]; yi_e=(]; 


7, Determine the user preferences and data file 
%Z names from the GUI. Go load the file from 
7%, the DATA directory. 
ipfnamei=get (loadi_ed,’string’); 
ipfname2=get (load2_ed, ’string’) ; 
edie... 
cd data 
eval([’load ’,ipfnamei,’.mat’]); 
cgrateA=cgrate; 
numfiles=1; 
if ( get(rb_file2,’value’) == 1 ), 
eval([’load ’,ipfname2,’.mat’]); 
cgrateB=cgrate; 
numfiles=numfiles+ti; 
end ; 
eae... 
cd model 


4 This is the target actuator rate: the level plane 
4 slicing the surface. 
udselect=str2num(get (limit_ed,’string’)); 


4 This records the user choice of either a 3D or 2D plot. 
if ( get(rb_3d,’value’) == 1) 
choice=1; 
else, 
choice=0; 
end; 


4 This captures the axis limits from the GUI. 
udmin=str2num(get (minrate_ed,’string’)); 
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udmax=str2num(get (maxrate_ed,’string’)); 
cgmin=str2num(get (mincg_ed, ’string’))/100; 
cgmax=str2num(get (maxcg_ed,’string’))/100; 
vhmin=str2num(get (nminvh_ed,’string’)); 
vhmax=str2num(get (maxvh_ed, ’string’)); 


4 This captures the view point preference. 
azim=str2num(get (azim_ed,’string’ )); 
elev=str2num (get (elev_ed,’string’)) ; 


_% The user may wish to compare to aircraft definitions 
4, side by side. Loop for 1 or 2 files. 
for comp=1:numfiles; 
if (comp==1), 
cgrate=cgrateA; 
else, 
cgrate=cgrateB; 
end; 


4, Preprocess the data by parsing by tail volume. 
[m,n]l=size(cgrate) ; 
num_vh=1; 
index (num_vh)=1; 
for 1=i:m-1, : 
if ( cprate@, De <Sepgratcas!, 1) 
index (num_vh+1)=i+t1; 
num_vh=num_vh+ti; 
end; 
end; 
index (num_vh+1)=m+1; 
file=’cgrate’ ; 
for 1=1i:num_vh, 
p=index(i) ; 
q=index(it1)-1; 
eval([file,int2str(i),” = "file, "@ qe eeape 
end; 


4 Curve fit the rate vs. cg data for each tail 

4 volume. The curve fit uses a simplex nonlinear 
4 fit. The general form of the function to fit 

4 to is two exponentials. 

global DATA 
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file2=’DATA’ ; 
tiles—'c’: 
fin={)> nlin=[]; 
for k=1:num_vh; 
eile? o— 9? file, int2str(k),’(:,2) ?,.... 
metewimucstr ik). ’(:,5)]; ’]) 


fone) t 0]; 
trace = Q; 
mol = .1; 


lambda = fmins(’myfit’,lam, [trace tol]); 

A = zeros(length(DATA(: ,1)),length(lambda)) ; 
for j = 1:length(lambda) 

A(:,j) = exp(-lambda(j)*DATA(:,1)); 


end 

Besaeeitales. = A\DATA(: ,2);?]): 
Mom=(lin cl; 

nlin=[nlin lambda] ; 
end 


4, Repeat the process, this time curve fitting 

4 for elev amp vs. cg for each tail volume. 
global DATA_E 

file2=’DATA_E?’ ; 

file3=’c’; 

lin_e=[(]; nlin_e=[]; 

for k=1:num_vh; 
Pumnotrabe2,? = [’ ,file,int2str(k),’(:,2) ’,.... 
imibemintostr k),?(:,6)): ’)) 


lam = [1 0]’; 
trace = Q; 
Gol =ac1; 


lambda = fmins(’myfit’,lam, [trace tol]); 

A = zeros(length(DATA_E(: ,1)) , length(lambda)) ; 
for j = 1:length(lambda) 

A(:,j) = exp(-lambda(j)*DATA_E(: ,1)); 

end 

epateGirires, = A\DATA_EC:,2);’1); 

ieinee line Cc]: 

nlin_e=[nlin_e lambda]; 
end 


4 Sample the functions determined above at uniform 
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4 values of cg. 
file4=’vh’; 
t=cgmin: .01:cgmax ; 
t=t’; 
for z=1:num_vh; 
for j = 1:length(t) 
r(j,z) = din(1,z)*exp(-nlinG@ez)+tG) ee 
lin(2,z)*exp(-nlin(2,z)*t(j)); 
e(j,z) = lin_e(1,z)*exp(-nlin_e(1,z)*t(j)) + 
lin_e(2,z)*exp(-nlin_e(2,z)*t(j)); 


end; 
eval([file4,’(:,zZ) = ° jfile int2ste)) (1. 1 eee 
ones(length(t),1); ’J) 

end; 


4 This is the resulting mesh. It could be plotted 

4, as is, but below we’1l decrease the mesh size 

/% using linear interpolation. 

ALLDATA=[] ; 

ALLDATA_E=[] ; 

for j=i:num_vh, 

ALLDATA=[ALLDATA t vh(:,j) r(:,j)]; 
ALLDATA_E=[ALLDATA_E t vh(:,3) e(:,3)]; 

end; 


7, Linearly interpolate the mesh between the values of 
4 tail volume. 
ht=[]; 
amp=[]J ; 
vol=[]; 
for k=1:((num_vh)-1); 
xi=vh(1,k):.01: (wh(1,k+1)-.01); 
x=vh(1,k:k+1); 
for j=1:length(t) 
y=r(j,k:k+1); 
y_e=e(j,k:k+1); 
yi(:,j)=interp1(x,y,xi); 
yi_e(:,j)=interp1(x,y_e,xi); 
end 
ht=[ht yi’]; 
amp=L[amp yi_e’]; 
vol=[vol;xi’]; 


end 

neta (: num vh) | ; 
amp=Lamp e(: ,num_vh)]; 
vol=[vol;vh(i,num_vh)]; 


7, Plot the data. 
if (comp==1), 
figure 
end; 
if (choice == 1), 
%, This is the 3D plot of the data. 
htt=ht’ ; 
[m,n]=size(htt) ; 
for p=1:m; 
hor Gg=1:n; 
Pemcnbibipe dg) 2 udmax | htt(p,q) < 0), 
ht (q,p)=nan; 
end; 
end; 
end ; 


if (comp==1), 

sunt (t*100,vol ,ht’) 

axis([cgmin*100 cgmax*100 vhmin vhmax udmin udmax] ) 
caxis([udmin udmax] ) 

grid; 

view(azim,elev) 


Jethis 1s thesrate jlamit slice. 
ipeSepmint O01: cemax-. 1; 

[mm ,nn]=size(ht) ; 
pln=ones(length(ts) ,nn) *udselect ; 
C=ones(length(ts) ,nn)*2; 
surface(ts*100,vol,pln’ ,C’); 
xlabel(’c.g. (%c)’) 
ylabel(’tail volume’ ) 
Zlabel(’actuator rate (deg/s)’) 
else, 

surf (t*100,vol,ht’) 
end; 


else, 


163 


7, This is the 2D intersection of the target rate limit 
4, and surface. 
vhsc=[]; 
cgsc=(]; 
[m,n]=size (ALLDATA) ; 
for j=1:3:n 
X1 = udselect; 
yi = spline(ALLDATA(: ,j+2) ,ALLDATA(: ,j) ,xi); 
vhsc=[vhsce ALLDATA(1,j+1)]; 
cgsc={cgsc yi]; 
end 


4 Fit a curve to the points representing the 2D graph. 
CGVH=[cgsc’ vhsc’]; 
[p, sl=polytit (COV 1) CGV Cyn) 


lam = [1 0]’; 
trace = Q; 
tol = .1; 


lambda = fmins(’myfit’,lam, {trace tol]); 

A = zeros(length(CGVH(: ,1)) ,length(lambda) ) ; 
for j = 1:length(lambda) 

A(C:,j) = exp(-lambda(j)*CGVH(: ,1)); 

end 

Eo= AVCGVHC. 2) 

lin_s=c; 

nlin_s=lambda; 


4 Sample the fitted function at uniform cg increments. 

vhsc=([]; 

cgsc=(]; 

t=cgmin: .05:cgmax ; 

t=t’; 

for j = 1:length(t) 

vhsc(j,1) = lin_s(1,1)*exp(-nlin_s(1,1)*t(j)) + 
lin_s(2, 1) *exp(—nlin_ s(2,1) +t). 

end; 

cgsc=t*100; 


4 Plot the 2D intersection of the level plane and surface. 
if (comp==1) , 
pleot( cose, vise, ou cose wise). 
eield:, 
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axis({cgmin*100 cgmax*100 vhmin vhmax]) ; 
remabel (7c... (4c) ’) 
ylabel(’tail volume’ ) 
else, 
plot(cgsc,vhsc,’*’ ,cgsc,vhsc) ; 
text (cgmax-.15,vhmint+.025,’* = 2nd File’); 
end; 

end; 


if (numfiles==2) , 

hold ; 

clear yi yl_e r e cgrate 
end; 


end; 


4 Return the user to the GUI directory. 
Cao 4 
cd gui 


4. THE DYNAMIC AIRCRAFT MODEL 
The model used by the Tail Sizing Design Tool simulates the longitudinal mo- 


tion of a conventional aircraft. The configuration may or may not include a canard. 
The upper shell of the dynamic model is either the Simulink file eom_b70r.m, for 
rigid-body dynamics, or the file eom_b70w.m, for flexible body dynamics. The trim 
and linearize routines supplied by MATLAB are use to determine the equilibrium 
point, and generate a linear model for the LMI problem formulation. The simulink 
shells above call the function eomb70r, or eomb70w, in order to compute the deriva- 
tive of the models’ states. The computation is based, in part, on the stability and 
control derivatives supplied by the user, which are retrieved via a nested function 
call explained later. The description of eomb70w follows, as the Simulink shell is self 
explanatory. 

Whhhhhhhhhhhhhhhhhhhhhhhhhhhhsh 


7, function: eomb70w.m i, 


hhhhhhhhhhhhhhhhhhhhhhhhhhhhhh 
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/, The vector passed to the function is the 
/ Vehicles ssbates: 
function accel = eomb70w(x) 


7, Retrieve the center of gravity , tail volume 
% and canard volume parameters 


load cgvh; 


/, Retrieve the normalized aircraft derivatives 
4 at the center of gravity and control volumes 
4 specified. Retrieve the flexible mode frequency. 


[rho,a,CDO,CL, CM, QF omegai mi, f1,s),b,e,m. tot. ita el| am 
b7Odataw(xcg,vh,split) ; 
load flexdata 


4 Parse the vector passed to the function into 

4 its components and use them to compute terms 

4 necessary to "denormalize" the force and moment 
/, coefficients. 

dc = x(1)5 de = x(2)j@dtrn— x@r 

u = x(4); w = x(5); q = x(6); 

theta = x(7), eta — x(6). etad — xo). 


/, Calculate the total velocity, vt, and dynamic 
4 pressure gravity (assumed constant) and 
4 angle of attack. 


Vite = (u°2 + w 2) 775. 
qbar = .5*rho*(vt*2); 
g =P 3272, 

alpha = asin(w/vt); 
normi = rhox*xvt~2*s/2; 
norm2 = rho*xvt*s*c/4; 
norm3 = rho*vt~2*s*c/2; 
norm4 = rho*evt*s*c*2/4; 
norm5 = rho*xvt~ 2*s*c/2; 
norm6 = rho*xvt*s*c*2/4; 


4 Based on the stability and control derivatives 
4 retrieved by the function {\em b70dataw.m}, the 
/, coefficients above and the vehicles’ states, 
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/, the total lift, drag, pitching moment and elastic 
4 generalized force on the vehicle is computed. 


L = normi*(CL(1) + CL(2)*alpha + CL(7)*de + CL(8)*dc + .. 
CL(5)*eta) + norm2*(CL(3)*0 + CL(4)*q + CL(6)*etad) ; 


Ln = L/(qbar*s) ; 
Dn = (CDO + 1.01*Ln~2/pi/a) ; 
D = Dn*qbar*s; 


M = norm3*(CM(1) + CM(2)*alpha + CM(7)*de + CM(8)¥*dc + .... 


CM(5)*eta) + norm4*(CM(3)*0 + CM(4)*q + CM(6)*etad) ; 


Qeta = norm5*(QF(1) + QF(2)*alpha + QF(6)¥*de + QF(7)*dc .... 


+ QF(4)*eta) + norm6*(QF(3)*q + QF(5)*etad) ; 


% The lift, drag and gravity force are resolved in 
4, the body reference frame. Then, the linear, angular 
% and elastc accelerations are computed in the body 


% reference frame. 


udot = -q*w - g*sin(theta) + (-cos(alpha)*D +.... 
sin(alpha)*L + To*dtrt)/m; 
wdot = q*u + g*cos(theta) + (sin(alpha)*D - .... 


cos(alpha)*L)/m; 
qdot = M/Iyy; 
etadd = -omegai(1)“2*eta - 2*omegai*Tixetad + Qeta/mi; 


Notice that the prior program utilized a nested function call to retrieve the 


total vehicle stability and control derivatives. These derivates are supplied by the 


function b70dataw, which provide user supplied values of the wing-body, tail and 


canard contributions to lift and pitch. If a flexible mode is retained, the user must 


also specify the mode shape, in vacuo frequency, and generalized mass of the first 


bending mode. Presented below is data for the generic HSCT model termed Ref A. 


When the flexible mode is included, the model is termed Ref B. 


elaltdstilatstcllaletate tlle late latrtetelataletet 
4, function: b70dataw.m i, 


167 


Ahhhhhhhhhhhhhhhhhhhhyhyhhhhhth 

74, To begin, the user must specify the physical 
4, description of the aircraft including the 

74, frequency and generalized mass associated 

%Z with the flexible mode. 


7, parameter/value/units/description 
32.174; % gravitational constant 
go = Cle, sipss2 

W = 480425; % lbs 

m= W/G; % slugs 

Iyy = 2.181e7; % slug-ft"2 

S = 0299; Apt 2 

C= 9c Oo eee 

cbar = ¢;%c 

b= 105; % wing span 

U= 422; % fps 


© 
i" 


rho = .0023769; % density 

Q = .S5*rho*U+U; % lof/ft*ft 

a = 1.751; % aspect ratio 

lh = 85; % tail moment arm to ac 
lc = 85; % canard moment arm to ac 
ih = 0; % tail incidence angle 

To = 1000; % nominal thrust lbs 
cloc = .05; % canard location 
tloc = .95; % aero center of tail location 
ve = split; 4% canard volume 

vh = vh; % tail volume 

omegai = 6.29; % flexible frequency 
mi = [500]; % generalized mass 

Ti = .02; &% flexible damping 


The stability and control derivatives of the 

canard, wing-body and tail are entered. Those listed 
below are for a generic HSCT "like" aircraft termed 
ref A. Following the development of Wykes’method, 
the user can enter terms that residualize the effects 
of the elastic motion. In this case, only the rigid 
body derivatives have been entered. The nomenclature 
uses either no extension or the extension (_r) to 
indicate that the term is a rigid body term. The 
extension (_e) is used for flexible body terms. 


a a a a 
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4 Finally, the extension (_eor) is used for terms that 
% scale the rigid body terms due to flexible effects. 


i Drag data 
CDO = 0.002; 
CDU = 0.0; 


yebatt data 

eljtt_ total = 0.161; 
clift_vh_term = 0.801642; 
clift_wb_afO_r = 0.04; 
clift_dwb_afO_er = 0.0; 
clift_wb_af_r = 0.0422; 
clift_wb_af_eor = 1.0; 
elitbowb_g_r = 0.0558; 


clift_wb_q_eor = 1.0; 
clift_wb_adf_r = 0.0; 
clift_wb_adf_eor = 0.0; 


clift_wb_n_e = 0.0; 


clift_wb_dq_e = 0.0; 
Qeicteh qr = 0.123; 
clift_h_gq_eor = 1.0; 
elatt oh adf_e = 0.0; 


clift_h_n_e = 0.0; 
Glitbon dg. e = 0.0; 
elattehn-ahO_e = 0.0; 
clift_h_ah_r = 0.0455; 
clift_h_ah_eor = 1.0; 
elatt el afO_r = 1.8; 
clift_e_afO_er = 0.0; 
clift_dah_dloadh = 0.0; 
Siitt ec afr = 0.71; 
clift_e_af_eor = 1.0; 


% Lift Data - Canard 
Gistteceqen «= 0.123; 
eliaite=ceq-eor = 1.0; 
clift_c_adf_e = 0.0; 
Glattec nue = 0.0; 
crintacedgie — 0.0; 
Slittec acl_e = 0.0% 
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clift_c ac _r “= 0.04555 


clift_c_ac_eor “= UE 
clift_dacsdloade = .0r 0. 
CLU = QO; 


/, Pitching moment data 
cm_total = -0.063618; 
EnO_wWwoer =] OcUr 
cmO_wb_er = 0.0; 
dcm_dcl_wb_r = 0.02; 
dcm_dcl_wb_er = 0.0; 
cm_wb_q_r = -0.0298; 
cm_wb_q_eor = 1.0; 
cm_wb_adf_r -0.0067; 
cm_wb_adf_eor = 1.0; 
cm_wb_n_e = -0.01; 
cm_wb_dq_e = 0.0; 

CMU = 0.0; 


4, Based on theses derivatives and the physical 
” characteristics of the aircraft, some useful 
4, coefficients are calculated. 


ki = cbhar/(2*U); 
k2 = cbar*vh/Lh;: 
k3 = S*cbar/lh; 
k2_c = cbhar*vc/lc; 
k3_c = S*cbar/lc; 


% Tail contribution coefficients 


Alpha_1 = clift_h_ahO_e... 
/(1-Cclift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
Alpha_2 = -(clift_h_ah_r*clift_h_ah_eor*.... 

(elift elat0ert+clittwcsatOmec) ae 
/(1-(clift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
Alpha_O = Alpha_1 + Alpha_2; 

Alpha_ih= clift_h_ah_r*clift_h_ah_eor... 
/(1-(Cclift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ) ; 
Beta = (clift_h_ah_r*clift_h_ah_eor*... 
(i-clift_e_af_r*clift_e_af_eor))... 
/(1-(clift_h_ah_r*clift_h_ah_eor*clift_dah_dloadh*Q*vh*k3) ); 
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7, Canard contribution coefficients. Note, the 
% canard is treated simply as a forward tail. 


mipma.oc- clift_c_acO_e... 

Pee (eclittec de rtcliftt_c_ac_eor*clift_dac_dloadc.... 
*Q*vc*#k3_c)); 

Adpha_ic= clift_c_ac_r+#clift_c_ac_eor... 

Pela clitt) c ac r#clift_c_ac_eor*clift_dac_dloadc.... 
*Q*Vc*k3_c)); 

Peeamem—mCclittec ac n*clift_c.ac_eor*(1-clift e_af r.... 
eelift_e_af_eor))... 

Pele (cliaft_¢_ac_r*clift_c_ac_eor*clift_dac_dloadc*Q*#vc*tk3_c)); 


7, Finally, the total vehicle stability and control 

7, derivatives are computed based on the canard, 

wing-body and tail contributions. The build up 

4% 1s based on a standard development presented in 

% many aircraft dynamics texts; here I used Etkin. 

4, The computation of the elastic stability and control 

74, derivatives follow the development earlier in this text. 


ss oS 


I 


Form static stability derivaties needed for trim. 


CEGe= clift_wb_af0O_r + clift_dwb_af0O_er + 
k2*Alpha_O + k2_c*Alpha_0Oc; 

CLA = 57.3*(clift_wb_af_r*clift_wb_af_eor + 
k2*Beta + k2_c*Beta_c); 

CLIS = k2*Alpha_ih; 

CLIC = k2_c*Alpha_ic; 

CMO = cmO_wb_r + cmO_wb_er + (dcm_dcl_wb_r + 
denmaciewo er) *(clitt_wb_afO_r + clift_dwb_afO_er).... 
- Alpha_0*vh +Alpha_Oc*vce + CLO*(xcg-.25) ; 

CMA = 57.3*((dcm_dcl_wb_r + 
wenmcelewomer) + (clitt wb at réclift wb af eor)... 
- Beta*tvh + Beta_c*vc) + CLA*(xcg-.25) ; 

CMIS = -Alpha_ih*vh + CLIS*(xcg-.285) ; 

CMIC = Alpha_ic*ve + CLIC*(xcg-.25); 


% Compute dynamic stability derivatives needed 
Pp y y 
Jeter linearization. 
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CLQ = 57.34#(clift_wb_q_r+ellittawpedecor tae 
clift_h_q_r*clift_hog_cor*tk2 -velit tee kqe rs eee 
cliftle_queon*kJec)- 

CLAD = 57.3*k1*(clift_wb_adf_r*clift_wb_adf_eor + 
clift_h_adf_e+*k2 + clift_eradt setk2ec)- 

CMQ = 57.3*((dcm_dcl_wb_r + dcm_dcl_wb_er)*.... 


(clift_wb_q_r*clift_wb_q_eor) + cm_wb_gq_r*cm_wb_q_eor -.... 


elift_h der elite oheqs@com aan 

clift_c_q_r*#clift_c_q_ eor+vc) 4CnGsGrece—e zoe 

CMAD = 57.3*((dcm_dcl_wb_r + dcm_dcl_wb_er)*.... 
(cm_wb_adf_r*cm_wb_adf_eor) 

- clift_h_adf_e*vh - clift_c_adf_e*vc) + CLAD*(xcg-. 25); 


4 Compute stability derivatives for the 

4 generalized elastic coordinates. When these are 
4 non-zero, the model is aeroelastic. The resulting 
4 model is termed Ref B is Chapter 2. 

[phi, phidot] =modeshpe( .25*250) ; 

CLeta = CLA*phidot; 

CLetad = CLA*phi/U; 

CMeta = CLA*(xcg-.25)*phidot ; 

CMetad = CLAD*phi*(xcg-.25)/U; 

QO = -CLO*phi; 

QA = -CLA*phi; 

QQ = -CLQ*phi; 

QETA = -CLA*phidot*phi; 

QETAD = -CLAD*phi*phi/U; 

[phi, phidot]=modeshpe(tloc*250) ; 

QIS = -CLIS*phi; 

[phi, phidot]=modeshpe(cloc*250) ; 

QIC = -CLIC*phi; 


7%, Combine the individual terms into a vector 
(eto be returned 

CL = [CLO-CLA- CLAD CLO Cletas¢letadecincec hl clr 
CM = [CMO CMA CMAD CMQ CMeta CMetad CMIS CMIC]; 
QF = [QO QA QQ QETA QETAD QIS QIC]; 


APPENDIX B. RFTPS SUPPORT 


1 INTRODUCTION 
This appendix documents the C-code unique to the RFTPS design that is not 


automatically generated via the autocode tools. The focus is on the serial module, 
which processes both the IMU and GPS data. More complete hardware descriptions 
of the I/O modules supported by the AC100 architecture can be found in the AC100 
Model C30 Supplemental Reference Manual. Wiring diagrams and pin assignments 


are on file in the Avionics Lab. 


2. SERIAL MODULE 


The serial module contains two multi-mode, bi-directional, channels. The user 
interface to the device is through the file user_ser.c. Within user_ser.c, three functions 
must be defined. The first, get SERIAL_parameters, is used to initialize the channels. 
The second, user SERIAL-out, defines the output from each channel and the third, 
user_SERIAL-_in, processes received data in the input buffer. Currently, the serial 
module is not used to send data. Some of the shell in user_ser.c is not shown. “The 


documentation focuses on RFTPS specific code associated with processing the GPS 


and IMU data. 


OOOO OG a I ak a a a ak ak ak ak ak ak 


**x File >: user_ser.c 

** Project RE DPS 

** Edit level : 4 

* 

**x Abstract: : File contains functions which the user 


** must define to interface with IP-SERIAL device 
** driver. 

** The functions must be called 

* 

** get_SERIAL_parameters 

** user_sample_SERIAL_in 

** user_SERIAL_out 
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kk 
xx Templates for the functions are provided. 

2k 

** get_SERIAL_parameters 

** Function sets the asynchronous communication 

*k parameters for the IP-SERIAL module. Ring buffer 

**k sizes used to store received data must also be 

xk specified. 

2k 

** user_SERIAL_out 

** Function 1s called every scheduler interval. The 

*k user 1s responsible for creating a byte stream from the 
*k models floating point outputs. The user must ensure 

xk that the when writing these bytes to the output buffers 
**x that the buffers are not overflown. 

2k 

xk user_sample_SERIAL_in 

xk Function is called every sampling interval. The 

xk user 1s responsible for filling the floating-point 

**k vector which is used as input to the model for 

xk the current sampling interval. 

2k 

2k 

** Modifications: 


kK -------------- $$ - $$ $= - = - - -- = -- = - = = == = 
eK Creation : 7-i--93 Henry Tominaga 

2k Revised : 8-23-93 Brent Roman 

2 K Revised : 11-18-93 Steve Lynch 

ok Revised : 9-1-94 Steve Lynch 

# Revised : 1-17-96 Eric Hallberg [New IMU] 

2K Revised : 3-15-96 Eric Hallberg [IMU Serial A 

* > / GPS Serial B] 

* > Revised : 5-01-96 Eric Hallberg [IMU binary] 

* / 


#include <stddef.h> 
#include <stdlib.h> 
#include "sa_types.h" 
#include "stdtypes.h" 
#include "actypes.h" 
#include "i1oerrs.h" 
#include "errcodes.h" 
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#include "1iodefs.h" 
#include “iodriver.h" 
#include "“ipserial.h" 
#include “printx.h" 


#define NULL 0 


/* semaphores and serial parameters for each 
physical channel */ 
private struct 
{ ISI_BOOLEAN allSent; ISI_BOOLEAN broken; 
unsigned baud;} line_state[2]; 


struct user_type 

{ 
int update_interval ; 
int update_count; 


i 


typedef struct _bytes 
{ 

unsigned bytel 
unsigned byte2 
unsigned byte3 : 
unsigned byte4 
}_bytes; 


we 


C © © @ 


typedef union float_char 


‘i 
leat eda; 
_bytes ch; 
ft oat echar- 


/* global variables used in user_ser_in ch A 
IMU data processing */ 


u_int buffer_data[40] ; 
int index; 

inte trek = 0: 

int first_frame_a = 0; 


175 


float lastetlodteaiia i. 
float last sfloate aan mane 
float last_float_a_2([14] ; 
float last_float_a_3[14]; 
tents missed_cr = 0; 


/* global variables used in user_ser_in ch B 
GPS processing */ 


int first_frame_b = 0; 
float last_float_b[50], sol=299792458.0, 11£0=1575420000.0, k2; 
unsigned long icpoldi, icpold2, icpold3, icpold4, 
icpoldS, icpold6; 
int num_bytes_old=0; 


/* Function: get_SERIAL_parameters ++++++++++++++4++++4+4+4 
** Returns: 
2 NONE 
a7 
public void get_SERIAL_parameters 
( 
unsigned int hardware_channel, 
volatile struct user_param *device_param, 
volatile struct ring_buffer_param *rec_buffer, 
IQdevice kdevice) 
{ 

struct user_type *user_ptr; 

serial_param_type *serptr; 


serptr = device->parameters; 


if (SERIAL_USER_PTR == NULL) 

{ 
SERIAL_USER_PTR = (void *)malloc(sizeof (struct user_type)) ; 
user_ptr = (struct user_type *)SERIAL_USER_PTR; 


user_ptr->update_interval = 1000; 
user_ptr->update_count = 0; 
} 


/* Both IMU and GPS sensors send their data at 
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9600 baud/RS-232 standard */ 


if (hardware_channel == chanA || hardware_channel == chanB) 
{ 

/* set parameters for channel A (IP-SERIAL channel 1)*/ 
device_param->parity = NONE; 
device_param->baud_rate = 9600; 
device_param->stop_bits = ONE; 
device_param->transmit_data_size = 8; 
device_param->receive_data_size = 8; 
device_param->clock_multiplier = 16; 


/* Set size for receive ring buffer. 

Here it is set to 10x the message size */ 
rec_buffer->buffer_size = 600; 

i 

else 

{ 

printx("INVALID CHANNEL\n") ; 

} 

} /* get_SERIAL_parameters */ 


/* Function: user_SERIAL_out +++++++++++++++4++4+4+4+4+4+4 
* / 
public RetCode user_SERIAL_out (IOdevice *device, 
tHloat model float, 
u_int ser_channel) 


i 
Byte ebiater | 201); 
u_int ah 
float_char data; 
RetCode return_val; 


return_val = OK; 

DCO OOO OO OO I kK Kk ak kK 

* Given floating point model output, 

* create buffer which contains bytes 

* to be transmitted across 

* serial channel. 

OOOO OK I a aK kak kak ak ak / 
if (numbytes_in_buffer(device->parameters) == 0) 
sf 


for (i = 0; i < 5; i++) 


late 


{ 
data.fl = model_float[il; 
cbuffer[0] =(Byte) data.ch.byte1; 
cbuffer[1] =(Byte) data.ch.byte2; 
cbuffer[2] =(Byte) data.ch.byte3; 
cbuffer[3] =(Byte) data.ch.byte4; 


[DOO OR OR RR OK IO a IK a 2 2k 

* Fill the output buffer with data to be transmitted 

* by the background portion of the serial driver 

FOO OO OO IC ok a a Rok a KK / 
return_val = write_serial (device->parameters ,cbuffer,4) ; 
if (return_val == -1) 
{ 
printx("Error: Buffer overflow in User_ser on 

output to serial\n") ; 


} 


} 


return OK; 
} /* user_SERIAL_out */ 


/* Function: user_sample_SERIAL_in +4++++++++++++++++ 
* / 
public RetCode user_sample_SERIAL_in(I0device *device, 
flea model_float{[], 
ene ser_channel) 
{ 
/*#* KA KKAA-KREADING CHANNEL A #2 0%% TMU ORK / 
if (ser_channel==chanA) 
{ 
/* Create a software buffer to hold receive 

buffer data */ 
lint sottabucrer loool 
/* Declare varibiables needed */ 

unsigned char item[1]; 

float k; 

Tic ey 20a 

int s = QO; 

int p = 0; 

u_int w_high; 

u_int w_low; 
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struct user_type *user_ptr; 


/* The IMU data is sent as twos complement 
and must be scaled*/ 
float scale[] = {180/8191.0, -2.0/8191.0, 
-2.0/8191.0, -2.0/8191.0, 
200.0/8191.0, 200.0/8191.0, 
200.0/8191.0, 180.0/8191, 
180.0/8191.0, 10.0/8191.0, 
10.0/8191, 10.0/8191.0, 
10.0/8191.0, 10.0/8191.0}, 
ieeaulicgedatall = 1 10.0, 10.0, -10.0, 10.0, 10.0, 10.0, 
MmMOmE Onn 10.0, 1.0, 1.0, 1:0, 1.0, 1.0}: 


DOO ORO OR RI A I aK aK aK a ok 3k 2k ak 2k 2k ok 

set user pointer to buffer allocated 

by get parameters 

this buffer is passed around with the 

structure device 

and should only be accessed via the SERIAL_USER_PTR 


x* + & © HF 


OOOO I RK I aK I ok ak aK 2k 2k 3 32k 2k 2k / 
user_ptr = SERIAL_USER_PTR; 


/* The first time the serial port is read; 
load default data */ 

if (first_frame_a==0){ 

for (j=0;j<14;j++){ 

last_float_a[j] = default_datalLj]; 

last_float_a_1[j] = default_datalj]; 

last_float_a_2[j] = default_datal[j]; 

last_float_a_3[j] = default_data[j]; 

}/* end if*/ 


model_float = last_float_a; 
index = 0; 

first_frame_a = 1; 

return OK; 

}/* end if*/ 


/* Process hardware buffer by moving data into 
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the software buffer */ 
while( (numbytes_in_buffer(device->parameters)) > 1 )f{ 
read_serial (device->parameters,1,item) ; 
soft_buffer[s]=(u_int) item[0]; 
s=st1; 
} /* end while */ 


/* The IMU data string ends with a carriage return. 
Always sync on the carriage return*/ 

while ( p < (s))¢f 

if (soft_buffer[p] == ’\r’){ 

missed_cr = 1; 

index = 0; 

}/* end if*/ 

else{ 

if (missed_cr == 1){ 

buffer_data[ index ] = soft_buffer|[p]; 

}/* end if */ 


index=index + 1; 
if (index == 30){ 
index = 0; 
missed_cr = 0; 
}/*end if*/ 

}/* end else*/ 
p=pti; 

} /* end while */ 


/* Implement some rudimentary error detection to remove 

* blatantly bad data from the IMU. Since the IMU does not 
send a check sum, use multiple buffers to compare present 
with past values. Remove data readings that are not 
physically possible. Send last good data. */ 


ee 


tick = tick + 1; 
if (tick == 4) 
tick = 1; 


if (tick == 1) 


{ 
for ( y=0; y<14; yt++){ 
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/* The number is in binary, 2’s complement. 
Convert to decimal and scale */ 
w_high = buffer_dataLy*2] ; 
w_low = buffer_dataLy*2+1]; 
k = 128*(w_high & Ox7F) + (w_low & Ox7F); 


me (k > 8191) 
k = k - 16384; 


last_float_a_ily] = k*scaleLly]; 


memcabs last ftloat_a_ily| - last_float_a_3[y]) < .1) 
last_float_aly] = last_float_a_ily]; 


} /* end for loop */ 
}/*end if*/ 


if (tick == 2) 
x 
for ( y=0; y<14; yt+){ 


w_high = buffer_dataLy*2] ; 
w.low = buffer_dataly*2+1]; 
k = 128*(w_high & Ox7F) + (w_low & Ox7F); 


erik > 8191) 
k = k - 16384; 


last_float_a_2ly] = k*scaleLy]; 


if (abs(last_float_a_2ly] - last_float_a_ily]) < .1) 
last_float_aLly] = last_float_a_2[y]; 


} /* end for loop */ 
}/*xend if*/ 


if (tick == 3) 
{ 
for ( y=0; y<14; y++){ 


w_high = buffer_data[y*2] ; 
w.low = buffer_dataLly*2+1] ; 


18] 


k = 128*(w_high & Ox7F) + (w_low & Ox7F); 


if (k > 8191) 
k = k - 16384; 


last_float_a_3[y] = k*scalel[y]; 


if (abs(last_float_a_3[yl) = lastefiloatean 2 ya ee) 
last_float_aly] = last_float_a_3ly]; 


} /* end for loop */ 
}/*end if*/ 


/* Pass the IMU data back up. 

for (j=0;j<14;jt++){ 
model_float{j] = last_float_a[j]; 
}/* end for*/ 

}/* end if ch a */ 


/***%*** READING CHANNEL B**** GPS 44K / 
if (ser_channel==chanB) { 


unsigned char itemb[1], msg1[600], checksumi[1], 
checksum2[1], mychecksum1 ,mychecksum2, status, 
mode1, mode2, mode3, mode4, mode5, modeé, 
msg3Ll45], checksum3[1], mychecksum3, msg4[62], 
checksum4[1], mychecksum4; 


SCrUCE USEFLtYpe =USeErLpLr: 


unsigned long 1, 3), 1eplsetep2eaaepsr 
icp4, icp5, icp6, 

ephi, eph2, eph3, eph4, 

ephS, ephé; 


long lat, lon, vel, heigps, heimsl; 


int numi, num2, num3, sati, sat2, 
sat3, sat4, sat5, sat6, 

satt1, satt2, satt3, satt4, 

satt5, satt6, num4, 


velnor, velup, veleast; 


float psrngi, psrng2, psrng3, psrng4, 
psrng5, psrng6, psrngratel, psrngrate2, 
psrngrate3, psrngrate4, psrngrated, 
psrngrate6, loctime, sattimel,sattime2, 
sattime3, sattime4, sattime5, sattime6, 
ecefx, ecefy, ecefz, velnori, veleasti, 
velupl, hd, timeref, pscor1l, pscor2, 
Pacors, pscor4, pscors, pscor6, 

Pareeor!, psreor2, psrcor3, psrcor4, 
Pencors, psrcor6, locsec, locsecf, 
satsecl, satsec2, satsec3, satsec4, 
satsec5, satsec6,satsecfl, satsecf2, 
satsecf3, satsecf4, satsecf5, 

satsecf6, heigpsi, heimsli, 

vel1, lati, loni, ecefx1l, ecefy1, ecefzl; 


DOO OO COG I OO aR I a a RK ak a kk ak ak 
* set user pointer to buffer allocated by get parameters * 
* this buffer 1s passed around with the structure device * 
* and should only be accessed via the SERIAL_USER_PTR * 
* define * 
OO Oc oo GI OK a I ka kak kk kk / 


Sere per = SERTAL_USER_PTR; 


[DRO OO ROR OR ORG RK RK Rk kK a Kk 
* set user pointer to buffer allocated by get parameters * 
* this buffer is passed around with the structure device * 
* and should only be accessed via the SERIAL_USER_PTR * 
* define * 
ROO COO oO a a aR ai a Ia aK aK a a / 


/* If first time read, load default data */ 

ie tipstetrame b=—0) + 

for (i=0;1<50;i++){ 

ie vetloateb ii |—7.0- 

first_frame_b=1; 

i1cpoldi=0.0; icpold2=0.0; icpold3=0.0; icpold4=0.0; 
i1cpold5=0.0; icpold6=0.0; 

}/*end for*/ 
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}/*end if*/ 


/* The message 1s 61 bytes long and comes in once a second. 
* Wait for the whole message to come in. */ 
numi = (numbytes_in_buffer(device->parameters) ) ; 
if ((numi > 60) & (numi1 == num_bytes_old) ){ 


for(i=0;i<numi ;it++){ 

read_serial(device->parameters,1,itemb) ; 

msg1(iJ=itembL0] ; 

}/* end for */ 

}/*end if message in*/ 

/* The message is in ASCII Hex and starts with @@BA. 
Search for start of string. */ 


for (i=4;i<num1;i++){ 


if (msgili-4]==’0’ & msgili-3]==’@’ & 
msgi[i-2]==’B’ & msgili-1]==’a’ ){ 


checksumil0] = msgi[it+é1]; 
/*Convert message bytes to decimal numbers*/ 


lat = msgi[i+11]*0x1000000 + msgi[i+12]*0x10000 + 
msgilit+13]*0x100 + msg1[i+14]; 


lon = msgili+15]*0x1000000 + msg1i[li+16]*0x10000 + 
msgili+17]*0x100 + msgili+18] ; 


heigps =msgili+19]*0x1000000 +msgi[i+20]*0x10000 + 
msg1lit+21]*0x100 + msgi[i+22] ; 


heimsl = msg1[i+23]*0x1000000 +msg1[i+24] *0x10000 + 
msg1[i+25]*0x100 +msg1i[i+26] ; 


vel = msgili+27]*0x100 + msg1i[i+28]; 


hd 


msg1li+29]*0x100 + msgi[30] ; 


hd 


hd 1070 
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status = msg1[i+60]; 

velli=vel/100.0; heigpsi=heigps/100.0; heimsli=heimsl/1i00.0; 
fati—1at/100.0; Voni=lon/100.0; 

/*GPS uses a check sum. Check if message is correct*/ 
mychecksum1=’B’~’a’ ; 

more GO }<01; j++) 4 

mychecksum1~=msg1i[litj] ; 


}/*xend for*/ 
/* If its correct, update GPS data. 
if (mychecksum1 == checksumi[0]){ 


Last_float_bLOJ=lati; 
tastat oat. bl1)]=ion1i: 
last_float_b[2]=heigps1; 
iastetiloat bil3)=heimsl1: 
fastetloat.bl4]l=vel1: 
Last_float_b[5]=hd; 
last_float_b([6]J=status; 

last _float_b[7]=checksum1 [0] ; 
last_float_b([8]=mychecksum1 ; 


}/*xend if checksum ok*/ 
1=1+60; 


}/* end if @@Ba */ 

/* Pass the GPS data back up */ 
for (j=0;j<8;j++)4 
model_float[j] = last_float_al[j]; 
}/* end for*/ 


}/* end if ch b */ 


return OK; 
} /* user_sample_SERIAL_in */ 


185 





[1] 


[2] 


(3 


el 


[4 


ee! 


[5 


[a 


(6 


Ns a! 


[7] 


[8] 


[9] 


[10] 


[11] 


[12] 


LIST OF REFERENCES 


I. Ashkenas, R. Magdaleno, D. McRuer, ”Flight Control and Analysis Methods 
for Studying Flying and Ride Qualities of Flexible Transport Aircraft,” NASA 
CR 172201, August, 1983. 


H. Ashley, Engineering Analysis of Flight Vehicles, Addison-Wesley, Reading, 
MA, 1974. 


Baxandall, P., and Liebeck, H., Vector Calculus, Clarendon Press, Oxford, 1986. 


Beaufrere, H., ” Limitations of Statically Unstable Aircraft due to the Effects of 
Sensor Noise, Turbulence, and Structural Dynamics,” Guidance, Navigation and 
Control Conference, Williamsburg, VA, August 18-20, 1986, pp. 721-731. AIAA 
PAPER 86-2203. 


Beaufrere, H., Soeder, 5., ” Longitudinal Control Requirements for Statically 
Unstable Aircraft,” NAECON 1986; Proceedings of the National Aerospace and 
Electronics Conference, Dayton, OH, May 19-23, 1986. Volume 2 (A87-16726 
05-01). New, York, Institute of Electrical and Electronics Engineers, 1986, pp. 
421-433. 


S. Boyd, L. El] Ghaoui, E. Feron, V. Balakrishnan, Linear Matriz Inequalities in 
System and Control Theory, SIAM, Phila., PA, 1994. 


Chilali, M., and Gahinet., P.,“ 7. Design with Pole Placement Constraints”, to 
appear in IEEE Transactions on Automatic Control, 1996. 


M. Chilali, P. Gahinet, C. Scherer, ” Multi-Objective Output-Feedback Control 
Via LMI Optimization,” IFAC Procedings, 1996, pp. 249-254. 


Costello, D., Kaminer, I., Carder, K., and Howard, R., “The Use of Unmanned 
Vehicle Systems for Coastal Ocean Surveys: Scenarios for Joint Underwater 
and Air Vehicle Missions,” Proceedings 1995 Workshop on Intelligent Control of 
Autonomous Vehicles, Lisbon, Portugal, March 1995. 


Craig, J. J., Robotics, Addison-Wesley, 1989. 


Douglas, I., and Koehler, K., “Low-Cost Aircraft Navigation System to Aid 
Global Climate Change Studies,” NASA News@mercury.hq.nasa.gov, RELEASE: 
96-99, May 14, 1996. 


Doyle, J., Glover, K., Khargonekar, P., and Francis, B., “State space solutions 
to standard Hy and H,, control problems,” [EEE Transactions on Automatic 


Control. Vol. AC- 34, No. 8, 1989, pp. 831-847. 


lon 


[13] Elgersma, M., Control of Nonlinear Systems Using Partial Dynamic Inversion, 
Ph.D. Thesis, University of Minnesota, Minneapolis, MN, 1988. 


[14] B. Etkin, L. Reid, Dynamics of Flight, Stability and Control, John Wiley and 
Sons, Inc., New York, NY, 1996. 


[15] Fryxell, D., Oliveira, P., Pascoal, A., Silvestre, C., and Kaminer, I., “Navigation; 
Guidance and Control of AUVs: An Application to the MARIUS Vehicle,” [FAC 
Journal of Control Engineering Practice, Vol. 4, No. 3, 1996, pp. 401-409. 


[16] P. Gahinet and A. Nemirovsku , “ General Purpose LMI Solvers with Bench- 
marks,” Proc. 32nd IEEE Conference on Decision and Control, pp. 3162-3165, 
San Antonio, TX, December 1993. 


[17] Gahinet, P., Nemirovski, A., Laub., A.J., Chilali, M., LMI Control Toolbox, The 
Mathworks Inc, May 1995. 


[18] Gill,P.E., Murray, W. and Wright, M.H., Practical Optimization, Academic Press 
Inc., N.Y -7N: Yo 19a pp ett s oles 


[19] Godhavn, J. M., “Nonlinear Tracking of Underactuated Surface Vessels”, 
preprint, 1996. 


[20] Healey, A., and Lienard, D., “Multivariable sliding mode control for autonomous 
diving and steering of unmanned underwater vehicles,” [EEE Journal of Oceanic 


Engineering, Vol. 18, No.3, pp. 327-339, 1993. 


[21] F.M. Hoblit, Gust Loads on Aircraft: Concepts and Applications, AIAA Edu- 
cation Series, American Institute of Aeronautics and Astronautics, Washington, 


DC, 1988, pp. 8-9. 


[22] R. Horn, C. Johnson, Matriz Analysis, Cambridge University Press, New York, 
NY, 1985. 


[23] Integrated Systems, Inc., “Xmath Interactive System Identification Module”, 
Part Number 000-0027-001, June 1994. 


[24] Kaminer, I., Motion Control of Rigid Bodies Using H. Synthesis and Related 
Theory, Ph.D. Thesis, University of Michigan, Ann Arbor, MI, 1992. 


[25] Kaminer, I., Hallberg, E., Pascoal, A., and Silvestre, C., “On the Design and Im- 
plementation of a Trajectory Tracking Controller for a Fixed Wing Unmanned 
Air Vehicle,” Proceedings 1995 American Control Conference, Seattle, June, 
1995. 


188 


[26] 


[27] 


[28] 


[29] 


[30] 


[31] 


[32] 


[33] 


[34] 


[35] 


[36] 


[37] 


I. Kaminer, R. Howard and C. Buttrill, "On the Development of Closed Loop 
Tail Sizing Criteria for an HSCT,” To appear in Journal of Guidance, Control 
and Dynamics. 


Kaminer, I., Pascoal, A., Khargonekar, P., and Coleman, E., “A Velocity Algo- 
rithm for the Implementation of Gain-Scheduled Controllers,” Automatica. Vol. 
3, 1995, pp. 1185-1191. 


J. Kay et al, “ Control Authority Issues in Aircraft Conceptual Design: Critical 
Conditions, Estimation Methodology, Spreadsheet Asessment, Trim and Bibli- 
ography,” VPI-Aero-200. 


Lin, C. Modern Navigation, Guidance, and Control Processing, Prentice-Hall, 
1991. 


E. Livne et al., ”Lessons from Application of Equivalent Plate Structural Mod- 
eling to an HSCT Wing,” Journal of Aircraft, Vol. 31, No. 4, July-Aug. 1994, 
pp. 953-960. 


A. Messac and P.D. Hattis, Physical Programming Design Optimization for 
High Speed Civil Transport,” Journal of Aircraft, Vol. 33, No. 2, March-April 
1996, pp. 446-449. 


A. Nemirovskii and P. Gahinet, “The Projective Method for Solving Linear Ma- 
trix Inequalities ,” Proc. 1994 American Control Conference, pp. 840-844, Bal- 
timore, MD, June 1994. 


R. J. Niewhoener and I. Kaminer, “On Integrated Aircraft/Controller Design 
Using Linear Matrix Inequalities,” Journal of Guidance, Control and Dynamics, 


Vol. 19, No. 2, March-April 1996, pp. 445-453. 


Papageorgiou, Evangelos, “Development of a Dynamic Model for a UAV”, Mas- 
ter’s Thesis, Naval Postgraduate School, March 1997. | 


Papoulias, F.,“ Stability Considerations of Guidance and Control Laws for Au- 
tonomous Underwater Vehicles in the Horizontal Plane,” Proceedings 7th In- 
ternational Symposium on Unmanned Untethered Vehicle Technology, Durham, 


N.H., 1991. 


Pascoal, A. M.,“ The AUV MARIUS: Mission Scenarios, Vehicle Design, Con- 
struction and Testing,” Proceedings 2nd Workshop on Mobile Robots for Subsea 
Environments, Monterey, California, May 1994. 


J.K. Ray, C.M. Carlin, and A.A. Lambregts, ” High-Speed Civil Transport Flight- 
and Propulsion-Control Technologies Issues,” NASA CR 186015, March 1992. 


189 


[38] 


[39] 


[40] 


[41] 


[42] 


[43] 


[44] 


[45] 


Roskam, J., Airplane Flight Dynamics and Automatic Flight Controls, Roskam 
Aviation and Engineering Corp., Ottawa, KS, 1979. 


Rugh, W. J., “Analytical framework for gain scheduling,” [EEE Control Systems 
Magazine, vol. 11, No. 1, 1991, pp. 74-84. 


Samson, C., “ Path Following and Time-Varying Feedback Stabilization of a 
Wheeled Mobile Robot,” Proceedings International Conference ICARCV’92, RO- 
[3 )) Singapore mio 


Schmidt, D.K, ”On the Integrated Control of Flexible Supersonic Transport Air- 
craft”, Guidance, Navigation and Control Conference, Baltimore, MD, August 


7-10, 1995, p. 258-265. AIAA-95-3200-CP. 


Silvestre, C., Pascoal, A. M., Fryxell, D., and Kaminer, I., “Design and Imple- 
mentation of a Trajectory Tracking Controller for an Autonomous Underwater 
Vehicle,” Proceedings 1995 American Control Conference, Seattle, June, 1995. 


Slattery, R. A., and Zhao, Y., “En-Route Descent Trajectory Synthesis for Air 
Trafic Control Automation,” Proceedings 1995 American Control Conference, 
Seattle, June 1995. 


Walsh, G., Tilbury, D., Sastry, S., Murray, R., and Laumond, J. P., “Stabilization 
of Trajectories for Systems with Nonholonomic Constraints,” JEEE Transactions 
on Automatic Control, Vol. 39, No 1, 1994, pp. 216-222. 


M. Waszak, D. Schmidt, ”Flight Dynamics of Aeroelastic Vehicles,” Journal of 
Aircraft, VOL. 25, NO. 6, June 1988, pp. 563-571. 


190 


INITIAL DISTRIBUTION LIST 


MG ence echnical Intormation Center .......caccc ccc ccc cvccccscccccee 
8725 John J. Kingman Road.,Ste 0944 
Ft. Belvoir, VA 22060-6218 


MPT mI TNNG@MO MO At, OF ca)s) « ceslelcinlg cs dic ss ces 0s 5 se cd ede elGine o cceeceouees 
Naval Postgraduate School 

411 Dyer Rd. 

Monterey, CA 93943-5101 


mebnemlcadc IiCamIner ssc. tcr es. sees eee 2s ee eee routes rs 
Department of Aeronautics and 

Astronautics, Code AA/KA 

Naval Postgraduate School 

Monterey, CA 93943-5101 


Mitts Mv EGTeM yi UU Gee = vc Gade cs eae ae eee a8 reeset 
Department of Aeronautics and 

Astronautics, Code AA/Ho 

Naval Postgraduate School 

Monterey, CA 93943-5101 


mr Eric N. Hallberg ...... seks ee eee 
10 Gillespie Ln 
Monterey, CA 93940 


191 





DUDLEY KNOX LIBRARY 


NAVAL POSTGRADUATE SCHOO. 
MONTEREY CA 93943-5101 





